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1.  SUMMARY 


The  purpose  of  the  program  of  research  reported  here  is  to  develop  and  demonstrate 
a  realistic  numerical  model  for  self-consistent  ionospheric  propagation  and  heating  of 
powerful  high  frequency  (HF)  wave';  by  conjoimng  existing,  separate  numerical  models  for 
high  frequency  wave  propagation  and  ionospheric  transport.  A  comprehensive  numerical 
implementation  of  plasma  geometric  optics  hns  been  enhanced  and  adapted  to  the  problem 
of  self-consistent  ionospheric  propagntion  of  high  power,  high  frequency  radar  signals. 
The  emphasis  in  this  work  is  on  the  reali‘?tic  and  rigorous  description  of  the  propagating 
waves  so  that  geometric  and  spectral  characteristics  of  nonlinear  ionospheric  effects  can 
be  predicted  in  detail.  The  plasma  geometric  optics  code  has  been  used  in  conjunction 
with  an  ionospheric  transport  code  to  determine  self-consistently  the  HF  radiation  field 
and  quasilinear  modifications  of  ionospheric  plasma  properties  due  to  ohmic  dissipation 
of  the  HF  radiation.  This  allows  detailed  study  of  nonlinear  phenomena  such  as  beam 
self-focusing,  formation  of  large-scale  density  cavities  and  striations,  and  induced  wave 
ducting.  Overviews  of  the  theory  underlying  the  existing  models  of  radiation  propagation 
and  ohmic  power  deposition  are  given,  and  enhancements  to  these  models  are  described. 
New  restilts  a^rhieved  using  the  enhanced  models  are  shown,  and  the  treatment  of  a  self- 
consistent  propagation  and  heating  problem  is  presented.  Ionospheric  propagation  effects 
on  HF  radar  signals  are  discussed,  and  possible  applications  of  the  self-consistent  modeling 
capability  to  radar  technical  issues  are  identified.  These  issues  are  especiaiUy  pertinent  to 
existing  and  planned  high-power  HF  radars. 

2.  INTRODUCTION 

The  ionosphere  is  an  inhomogeneotis,  amsotropic,  and  dissipative  medium  which  is 
frequently  used  as  a  propagation  channel  by  electromagnetic  waves.  When  the  wave 
power  is  low,  little  interaction  with  the  plasma  occurs,  and  the  propagation  is  quedified  as 
hnear.  When  the  wave  power  is  high,  the  wave  inter£u:ts  with  the  medium  and  changes 
its  properties,  and  a  nonlinear  theory  is  required  to  describe  the  process.  This  nonlinear 
action  can  be  either  disruptive  to  an  operating  system  using  the  ionosphere,  or  it  can  be 
intentionally  produced,  as  is  done  in  ionospheric  modification  experiments. 

In  the  HF  frequency  regime,  nonlinear  phenomena  have  been  identified  both 
theoretically  and  experimentally.  They  can  be  categorized  according  to  their  length 
2Lnd  time  scales  into  two  distinct  regimes:  large  scale  electromagnetic  phenomena, 
occurring  on  scales  on  the  order  of  kilometers  and  seconds,  and  small  scale  electrostatic 
phenomena,  occurring  on  scales  of  the  order  of  centimeters  and  milliseconds.  It  is 
useful  to  assess  the  contribution  of  each  phenomenon  to  the  nonlinear  HF  heating  of 
the  ionosphere  and  to  determine  which  mechanism  is  more  likely  to  dominate  imder  given 
HF  parameters  and  ionospheric  conditions.  For  example,  parametric  instabilities  have 
long  been  thought  to  play  a  major  role  in  plasma  heating^^l.  Specifically,  two  major 
instabilities,  parametric  decay  and  thermal  self-focusing,  have  been  frequently  invoked  to 
account  for  experimental  results.  Although  the  former  belongs  to  the  category  of  small 
scale  electrostatic  phenomena^^^ ,  the  latter  is  observed  as  large  scale  structures  (1-5 
kilometers)  drifting  across  the  ionosphere  with  a  time  scale  of  the  order  of  one  minute^^l . 
Different  theories  of  thermal  self-focusing  provide  different  estimates  of  scale  length  for  the 
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resulting  striations,  varying  from  meters  to  kilometers!^}.  To  date,  the  relative  contribution 
of  electrostatic  and  electromagnetic  phenomena  in  producing  the  observed  striations  is  still 
unclear.  More  recently,  strong  turbulence  theory  of  caviton  collapse  and  burnout!^}  has 
been  used  with  new  experimental  data  to  put  forward  a  mechanism  for  early-time  HF 
energy  absorption  which  is  finding  favor  over  the  more  traditional  assumption  of  weak 
turbulence.  Because  of  the  coUisional  nature  of  the  ionosphere,  it  is  estimated  that  a 
large  amount  of  the  energy  carried  by  the  HF  waves  is  dissipated  ohmically!^},  first  in 
the  highly  coUisional  D  and  E  regions,  and  then  in  the  F  region.  In  the  F  region  ohmic 
heating  produces  plasma  transport  primarily  along  the  magnetic  field  lines.  This  can 
create,  among  other  things,  large  scale,  field-aligned  thermal  density  cavities  which  can 
also  cause  self-focusing  and  can  induce  wave  ducting. 

For  aU  the  processes  mentioned  above,  a  self-consistent  analysis  which  relates  the  HF 
waves  emanating  from  the  transmitter  antenna  to  the  local  region  in  the  plasma  where  the 
energy  is  deposited  is  crucial.  Such  an  analysis  permits  a  true  estimation  of  the  pump  field 
amphtude,  polarization,  and  phase  needed  to  initiate  the  instabilities  discussed.  A  self- 
consistent  analysis  also  traces  the  different  mechanisms  to  a  common  HF  beam  injected 
into  the  ionosphere  so  that  the  relative  importance  of  the  competing  processes  can  be 
correctly  assessed.  Calculation  of  wave-plasma  interaction  processes  can  be  made  more 
meaningful  in  the  context  of  consistent  wave  and  background  plasma  information.  The 
results  provide  valuable  information  for  design  and  performance  evaluation  of  HF  radars. 

In  the  current  work,  these  issues  have  been  addressed  by  enhancing  and  adapting 
a  numerical  implementation  of  plasma  geometric  optics,  initiaUy  developed  at  Yede 
University!^},  which  uses  the  geometric  optics  approximation  to  calculate  wave  amplitude, 
phase,  and  polarization  in  inhomogeneous,  anisotropic,  three-dimensional  media.  This  code 
is  applied  to  the  specific  problem  of  oblique  ohmic  heating  of  the  ionosphere  by  powerful 
HF  waves  such  as  those  generated  by  HF  radar  transmitters.  The  power  dissipation  due  to 
the  geometric  optics  wave  field  is  used  to  determine  the  heat  soTirce  in  a  plasma  trEuisport 
code  which  predicts  ohmic  heating  of  the  ionosphere.  This  heating  in  turn  affects  the 
dielectric  tensor  of  the  ionosphere  and  the  propagation  of  waves  is  influenced.  A  block 
diagram  outlining  the  procedure  for  self-consistent  analysis  of  ionospheric  HF  heating  is 
shown  in  Figure  2-1. 

The  plasma  dielectric  tensor  is  implemented  as  a  separate  module  invoked  by  the 
geometric  optics  code.  This  modular  aspect  permits  the  use  of  various  types  of  backgroimd 
plasma  models  in  the  radiation  calculations.  For  the  current  effort  an  inhomogeneous, 
magnetized  and  dissipative  cold  plasma  model  is  used  for  the  dielectric  tensor.  The 
necessary  ionospheric  parameters  are  obtained  from  existing  standard  models. 

In  Section  3,  a  description  of  the  theory  underl^dng  the  plasma  geometric  optics 
code  is  presented.  Section  4  indicates  improvements  which  have  been  made  in  the  wave 
propagation  model,  including  the  models  used  for  the  background  ionosphere.  Some 
results  of  computations  for  oblique  ionospheric  propagation  in  the  presence  of  realistic 
ionospheric  parameters  are  presented  in  Section  5.  Section  6  gives  an  overview  of  the 
existing  ionospheric  transport  model.  In  Section  7,  upgrades  to  the  ionospheric  transport 
code  which  extend  its  applicability  into  the  lower  ionosphere  are  explained.  In  Section 
8  the  procedures  used  to  couple  the  wave  propagation  and  transport  codes  are  set  forth. 
Both  linear  and  self-consistent,  nonlinear  heating  calculations  are  presented  in  Section  9. 
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Section  10  offers  a  summary  of  propagation  effects  on  radar  sign^.  Overall  conclusions 
are  ^ven  in  Section  11,  along  with  anticipation  of  specific  applications  of  this  work  to 
radar  performance. 

Figure  2-1.  Block  Diagram  for  the  HF  Ionospheric  Propagation  Code 
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3.  OVERVIEW  OF  THE  RADIATION  DESCRIPTION 

Plasma  geometric  optics  is  an  asymptotic  theory  of  linear  electromagnetic  waves 
in  plasmas  with  weak  space  stnd  time  dependence  based  upon  the  well-known  WKBJ 
approximation  techniques.  Despite  the  relative  simplicity  of  the  computational  prescription 
resulting  from  this  as3rmptotic  theory,  complications  associated  with  anisotropy  of  plasma 
media  and  vectorial  wave  fields  reqmre  that  nontrivial  li('lds  be  generated  by  nmnerical 
rather  thzin  analytical  means.  A  technique  known  as  rny  tracing  is  used  for  numerical 
description  of  high  frequency  wave  propagation  in  the  geometric  optics  limit. 

Although  many  ray  tracing  computations  have  been  performed  and  reported^®!,  the 
code  utilized  here  is  the  first  to  determine  fully  the  variation  of  the  amplitude  and 
polarization  of  the  field.  The  code  generates  boundary  data  describing  transmitted 
radiation  in  a  form  compatible  with  the  asymptotic  theory  of  propagation,  launches 
geometric  optics  rays,  decomposes  them  into  geometric  optics  modes,  and  propagates 
each  of  these  along  separate  refracted  ray  trajectories.  The  amplitudes  and  p>olarizations 
of  the  radiation  modes  are  calculated  along  the  way,  which  is  essential  for  the  correct 
determination  of  power  deposition  resulting  from  coherent  superposition  of  several  waves. 
The  code  is  compatible  with  an  arbitrary  model  of  inhomogeneous,  anisotropic  background 
plasma  properties  as  they  contribute  to  the  local  plasma  conductivity  tensor,  since  roots 
of  the  local  plasma  wave  dispersion  relation  are  tracked  through  parameter  space  by 
performing  computations  directly  with  the  components  of  the  conductivity  tensor  itself. 
In  addition,  the  code  admits  general  weak  space  and  time  variation  of  the  background 
plasma  and  fields.  This  latter  property  allows  the  direct  prediction  of  spurious  radar 
doppler  shifts  caused  by  time-dependent  density  fiuctuations  in  the  ionosphere.  The  code 
detects  the  occurrence  of  caustics  and  implements  the  asymptotic  matching  which  results 
from  a  local  boundary  layer  analysis,  permitting  graceful  continuation  of  the  geometric 
optics  calciilation  as  shown  in  Figure  3-1. 

For  ionospheric  heating  calciilations,  this  capability  must  be  extended  beyond  the 
limitations  of  geometric  optics  so  that  accurate  field  values  are  obtained  everywhere  near 
the  caustic  surfaces.  This  may  be  accomphshed  in  a  straightforward  manner  a  posierioTi 
by  patching  the  geometric  optics  amphtude  in  the  neighborhood  of  the  caustic  singularity 
with  an  amplitude  obtained  from  an  appropriate  boxmdary  layer  analysis.  Birefringence 
and  Faraday  rotation  are  treated  in  a  fully  automated  and  general  fashion.  An  appropriate 
modification  of  the  asymptotic  theory  is  automatically  implemented  in  nearly  degenerate 
environments,  where  two  different  modes  have  nearly  the  same  local  dispersion  relation. 
Figure  3-2  illustrates  the  results  of  this  procedure  when  an  electromagnetic  wave  enters  a 
magnetized  plasma. 

This  procedure  results  in  a  correct  accounting  for  linear  mode  conversion  in  weakly 
anisotropic  plasmas.  It  zdso  allows  the  prediction  of  large  scale  spatial  modulation  of  wave 
dissipation  rates  due  to  long  wavelength  Faraday  rotation  in  the  lower  ionosphere. 

The  asymptotic  theory  upon  which  the  geometric  optics  code  is  based  may  be 
sxunmarized  as  follows.  An  equation  governing  linear  waves  in  a  magnetized,  weedcly 
nonuniform  plasma  follows  directly  from  Maxwell’s  equations  and  an  appropriate  linear 
constitutive  relation. 
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K  the  dependence  of  £  on  its  last  two  arguments  is  weak  then  nearly  plane  wave  solutions 
to  (1)  may  be  sought  of  the  form 

£(r,«)  =  £o(r.t)e''^<^''>  (2) 

where  the  complex  amplitude  vector  is  slowly  varying  in  space  and  time  and  the  real 
phase  function  tjfiaa,  nearly  linear  function  of  space  and  time.  When  the  aruatz  (2)  is  used 
in  the  wave  equation  (1)  there  results 

*A  •  £o  =  i(£o)  (3) 

accurate  through  first  order  in  small  inverse  space  and  time  scales.  This  approximation 
neglects  effects  which  are  of  second  order  in  the  ratio  of  wavelength  to  the  length  scale  of 
spatial  inhomogeneity  and  also  those  which  are  of  second  order  in  the  ratio  of  wave  period 
to  the  time  scale  of  temporal  nonuniformity.  Here  the  operator  L  may  be  written 

i(£o)  =  ■  £o)  + 


-iv-(aiA^.£o)-5^-(4A^-£o)-AA-£o  (4) 

where 

k  s  uj  =  (5) 

The  operators  with  hats  are  tmderstood  to  act  only  on  explicit  space  and  time  dependences, 
whereas  those  without  hats  act  on  the  implicit  dependences  through  k  and  a;  as  well.  If 
the  plasma  is  perfectly  uniform  then  the  right  hand  side  of  (3)  vanishes.  The  tensor  A  has 
the  form 

A  =  (^  -  (6) 

where  £  is  the  uniform  plasma  conductivity  tensor  evaluated  using  the  local  space- 
dependent  and  time-dependent  plasma  properties.  For  HF  waves,  ions  contribute  negligibly 
to  the  dielectric  properties  of  the  ionosphere  and  thermal  dielectric  effects  are  also 
ummportant.  Thus,  the  ionosphere  is  modeled  as  a  cold,  magnetized  electron  plasma 
with  Krook  type  dissipation  governed  by  collision  frequency  i/c,  for  which  the  conductivity 
tensor  has  the  form 


where 


(I)  =  u;  -I-  iz/c 


,  2^  =  — s,  0,2  = 

me  ^ 


47re^n 


and  R  is  the  backgroimd  magnetic  field,  e  is  the  magnitude  of  the  electron  charge,  m  is 
the  electron  mass,  c  is  the  speed  of  light,  and  n  is  the  electron  deiisity.  The  tensor  AA 
vanishes  in  the  cold  plasma  limit,  and  in  general  is  not  simply  expressible  in  terms  of  the 
tensor 
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A  series  solution  to  (3)  is  developed  based  upon  the  smallness  of  the  right  hand  side, 
whence 

>A-4°’=0  (9) 

and 

iA'4*’=i(4'”)-  (10) 

In  order  for  the  phase  function  tj;  to  be  real  when  the  conductivity  tensor  £  is  not  completely 
antihermitian,  an  additional  small  term  may  be  consistently  introduced  in  both  zero  and 
first  ordo:  which  contains  a  free  parameter  v,  so  that  the  system  (9)-(10)  becomes 

t'A  •  =  -iud^A  •  (11) 

and 

■A  •  4'* = i(4°’) + ■  4°’-  (12) 

The  zero  order  eqiiation  (11)  has  nontrivial  solutions  only  when 

dei(A  +  it'duA)  =  0  (13) 

which  defines  a  local  dispersion  relation 

w  =  n(fc,r,t),  (14) 

a  dissipation  parameter  v{k,r,t),  and  a  polarization  vector  H  such  that 

(A  +  tt'^A)  •  iZ  =  0  (15) 

where 

4°’  =  (16) 

and  Q  is  a  complex  scalar  amplitude.  The  polarization  vector  II  has  an  adjoint  V_  such 
that 

+  =  (17) 

When  the  first  order  equation  (12)  is  projected  onto  Y.  there  results,  after  some 
manipulations, 

^o  +  ’  Va  +  -(V  •  5^11)0  =  (18) 

where 

^  s ./+  •  (va.A)  -  i(a,aL,A) 

+\[c^u.f  ■■  OiA^)  •  r  -  ( w)' :  OiA)  •  m 

—miLf  ■■  (oLA^)  ■  £*  -  (diZ)^  ■  (0l4)  ■  IL] 

+i[£t .  .  £7  -  i^t .  .  Q 

-i:'.(AA).£r}  (19) 
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and  A  is  that  singular  value  of  ^+ti/5L)A  which  vanishes  when  (13)  is  satisfied.  Application 
of  the  method  of  characteristics  to  (18)  yields  the  ordinary  differential  equation 

5ra  +  i(V  .  d^Sl)a  =  (20) 

with  characteristic  curves  defined  by 

arr  =  d^.  (21) 

Here  r  is  a  time  variable  which  parameterizes  the  characteristic  curves.  Consequently  its 
value  is  propagated  through  space  by  the  characteristic  curves  and  it  may  be  considered 
as  a  curvilinear  space  coordinate.  Differentiation  of  (13)  leads  to 

=  _vn  (22) 

which  formally  completes  the  system  of  ray  and  amplitude  equations.  The  polarization  of 
the  wave  is  trsmsported  along  the  ray  by  differentiating  ( 15)  with  respect  to  r  and  solving 
for  the  derivative  of  12,.  Finally,  if  coordinates  and  S2  are  defined  which  parsimeterize 
the  boundary  surface  from  which  rays  are  laimched,  then  it  c£m  be  shown  that 


drj  =  (V  •  dkQ)j 

(23) 

where 

y  S  5*,  r  X  da^r.  ‘  drL. 

(24) 

Equation  (23)  can  be  used  to  simplify  (20),  which  becomes 

(25) 

where 

=  :  (8iA^)  ■  r  -  (VZ)'  :  OiA)  ■  a 

■■  ■  r  -  (ftz)' :  OwA)  ■  m 

+l\y}  (^a,A)  ■  m 

-i:''(aa)e:}.  (26) 

The  divergence  of  the  group  velocity  appearing  in  (23)  becomes  infinite  on  a  caustic  surface. 
In  practice,  therefore,  (25)  is  integrated  along  with  the  ray  equations,  the  Jacobian  J  is 
calculated  from  (24),  and  equations  determining  the  partial  derivatives  in  (24)  are  also 
integrated, 

drE  =  E-  +  R  •  i^kdk^)  (27) 

and 

= -£  •  ( VVD)  -  2  (dt  VQ)  (28) 
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where  ^  is  any  one  of  the  partial  derivative  vectors  appearing  in  (24)  and  2  is  the 
corresponding  partial  derivative  of  k-  The  advantage  of  using  (25),  (27),  and  (28)  is 
that  all  of  the  coefficients  and  dependent  variables  are  bounded  on  a  caustic  surface. 

When  the  background  is  isotropic  01  nearly  isotropic  the  foregoing  analysis  breaks 
down  and  the  analysis  must  be  repeated  for  the  case  in  which  two  singular  values  of  the 
wave  equation  tensor  vanish  or  are  both  small.  One  result  of  switching  between  the  two 
formalisms  is  the  bifurcation  of  rays  on  entry  into  nn  anisotropic  medium. 

The  ray  tracing  code  described  by  the  block  diagi  am  in  Figme  3-3  has  been  constructed 
using  the  formalism  described  above.  Examples  of  its  capabilities  are  displayed  in  Figure 
3-4. 

Figure  3-1.  Handling  of  a  Caustic  Surface  by  the  Plasma 
Geometric  Optics  Code 
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Figure  3-2.  Treatment  of  Birefiringence  by  the  Plasma  Geometric  Optics  Code 
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Figure  3-3.  Block  Diagram  for  the  Plasma  Geometric  Optics  Code 
and  the  Transport  Code  Showing  Main  Modules^  Input  and  Output 
Routines,  and  Counline  Interfaces 
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Figure  3-4.  Examples  of  Plasma  Geometric  Optics  Code  Capabilities 
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4.  ADAPTATION  OF  THE  WAVE  PROPAGATION  COMPUTER  CODE 
TO  IONOSPHERIC  MEDIA 


4-1.  Antenna  Pattern  and  Radiation  Launching 

Electromagnetic  fields  at  the  botmdary  surface  from  which  rays  are  latmched  must 
conform  to  the  assumptions  of  slowly  varying  amplitude  and  nearly  linear  phase  which 
tmderlie  the  geometric  optics  theory.  Fields  radiated  from  a  transmitting  antenna  mtist 
therefore  be  expressed  in  their  far  field  form  on  the  boundary  surface.  The  far  field  form 
in  conjimction  with  the  local  dispersion  relation  serves  to  determine  the  geometric  optics 
phase,  amplitude,  polarization,  and  associated  properties  of  the  transmitted  wave  on  the 
exterior  of  a  sphere  surrounding  the  antenna. 

An  idealized  concept  of  an  HF  radar  antenna  has  been  used  to  develop  a  simple, 
approximate  description  of  a  far  field  antenna  radiation  pattern.  This  description  has  been 
employed  to  determine  initial  data  for  the  ray  tracing  calculations  of  HF  propagation.  The 
antenna  is  conceived  of  as  a  rectangle  on  which  a  uniform,  monochromatic  surface  current 
has  been  excited.  The  tiniformity  of  the  surface  current  may  be  thought  of  as  the  result 
of  a  Kirchoff  type  of  approximation.  The  far  field  form,  derived  below,  is  used  to  compute 
amplitude  and  polarization  boxmdai^'  data  for  a  geometric  optics  field  calculation  on  the 
exterior  of  a  sphere  surrotmding  the  antenna.  The  sphere  is  chosen  to  be  large  compared  to 
the  antenna  but  does  not  extend  into  the  ionosphere.  Geometric  optics  rays  are  launched 
from  the  surface  of  this  sphere  and  the  prescribed  amplitude  and  polarization  data  are 
transported  into  the  ionosphere.  The  arrangement  is  shown  in  Figure  4-1. 


The  magnetic  field  associated  with  a  surface  current  K.  having  harmonic  time 
dependence  is 


(1) 

where  G  is  the  Green’s  function 

f.ikR 

(2) 

and  k  =  ujc.  In  the  far  field, 

e’*’’  1.  j 

VG  -  ik - 

r 

(3) 

where 

k  = 

In  this  case  the  magnetic  field  due  to  a  uniform  current  K  becomes 

(4) 

£(r)  =  — i  X  K  L'^l' 

CT  JS 

(5) 

This  integral  may  be  immediately  evaluated  over  a  rectangular  antenna  siirface  with 
dimensions  Ax  and  Ay  centered  at  n'  =  0  to  give 


&{l)  =  AxAy - 1  X  K 

CT 


j»m(t*xAx/2)  stn(l;yAy/2) 
)t,Ax/2  JkyAx/2 


(6) 
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where  kx  =  kx/r  and  ky  3=  ky/r.  The  corresp<mding  electric  field  is  given  by 


£=--4x£ 

u> 


=  -AxAy'^k  X  k  X 

cjr  kxAx/2  kyAx/2 


(7) 


These  fields  have  a  main  lobe  and  sidelobe  patterns  in  azimuth  and  elevation.  The  average 
power  emanating  from  the  antenna  is 


p2 

P  =  c-^AxAy 
8ir 


(8) 


where  Qq  is  the  magnetic  field 
current  by 


whence 


Thus  the  radiated  power  is 


at  the  antenna  surface.  This  field  is  related  to  the  surface 


— K  *  n  X  £0 

c 

(9) 

Bl  = 

c 

(10) 

P  =  —  AxAyif^ 
c 

(11) 

and  the  magnitude  of  the  surface  current  K  may  be  chosen  to  be  consistent  with  a  certain 
radiated  power  when  the  radiated  electric  fidd  is  evaluated  on  the  spherical  boundary. 


4-2.  Evaluation  of  Wave  Power  Deposition 

The  procedure  for  evaluating  the  radiation  electric  field  and  associated  power  flux  and 
volumetric  power  deposition  once  the  ray  tracing  results  are  in  hand  is  not  trivial.  E)ach 
geometric  optics  ray  mtist  be  considered  as  a  member  of  a  family  of  rays  which  describes 
a  single  geometric  optics  mode  j  having  the  form 

£j=£(i,(r,  (12) 

where  the  amplitude  vector  is  nearly  constant  and  the  phase  function  is  nearly  linear. 
The  transport  of  the  wave  phase,  amplitude,  and  polarization  along  a  ray  constitutes 
evaluation  of  a  single  geometric  optics  mode  along  the  ray  path.  In  numerical  practice, 
values  of  these  wave  parameters  are  available  at  discrete  points  along  the  ray  path.  Thus, 
for  each  geometric  optics  mode  there  is  a  three  dimensional  grid  of  points  at  which  the 
wave  parameters  are  available.  This  grid  is  generated  by  tracing  rays  laimched  from  a  two 
dimensional  boimdary  surface.  The  field  of  each  mode  at  an  arbitrary  point  in  space  and 
time  must  be  determined  by  a  three  dimensional  interpolation  procedure.  The  procedure 
is  as  follows.  Let  si  and  S2  be  coordinates  parameterizing  the  boundary  surface,  and  let 
r  be  a  coordinate  parameterizing  each  ray.  Then  wave  parameters  for  a  given  mode  may 
be  made  available  on  a  regular  grid  in  the  three  dimensional  ray  space  (si,  S2,  r)  after 
rays  have  been  traced.  Likewise,  the  positi<m  coordinates  of  these  gridpoints  are  tabulated 
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^  Suppose  that  the  field  of  that  mode  is  desired  at  a  point  £o-  ^h'st  it  is 

necessary  to  determine  ray  space  coordinate  values  (siq*  ^20^  '^)  corresponding  to  the  field 
point  £0  such  that 

no  =  r(«io»«20iTo)-  (13) 

This  mxist  be  done  by  inverse  interpolation,  interpolation  is  justifiable,  a  straightforward 
three  dimensional  numerical  rootfinding  scheme  will  solve  this  problem  with  any  reasonable 
initial  guess.  In  particular,  a  Newton  iteration  based  on  the  relation 

6r  =  6s\da^z  +  Ss2daiL  +  SrdrL  (14) 

where 

^£  =  2^)--r(5l,S2>^)  <15) 

works  well.  Here  the  partial  derivatives  can  be  evaluated  using  tabulated  ray  coordinate 
data.  The  increment  6t,  for  example,  may  be  determined  by  taking  the  iimer  product  of 
the  total  variation  6r  above  with  the  vector  dg^r  x  This  vector  is  perpendictdar  to 
both  ds.L  and  and  the  result  is 

6t  =  {^a^r  X  da^r)  •  6rl{daiT  x  d^jr)  •  drL.  (16) 

Increments  and  may  be  similarly  determined.  The  coordinates  s^,  S2,  and  r  are 
repeatedly  incremented  imtil  they  reach  sio,  S20)  und  tq. 

Having  determined  s^q,  S20?  and  tq  in  this  fashion,  the  wave  parameters  of  the 
geometric  optics  mode  may  be  forward  interpolated  in  three  dimensions.  That  is,  the 
phase,  amplitude,  and  polarization  at  sio«  ^20*  and  may  be  found  from  corresponding 
values  at  nearby  points  in  the  three  dimensional  ray  space  grid  determined  by  the  ray 
tracing.  The  entire  process  of  inverse  interpolation  followed  by  forward  interpolation  must 
be  repeated  for  each  geometric  optics  mode  which  passes  across  the  field  point.  In  a 
typical  oblique  propagation  scenario  there  may  be  four  modes  passing  through  a  field 
point:  upward  laimched  and  downward  refiected  instances  of  ordinary  and  extraordinary 
magnetized  plasma  waves.  Coherent  superposition  of  these  waves  is  necessary  in  order  to 
display  such  featmres  as  Faraday  rotation  of  the  polarization  vector  and  associated  spatial 
oscillation  of  the  power  deposition. 

Time-averaged  conservation  of  energy  for  complex,  harmonic  fields  is  expressed  as 

V-5=-p  (17) 

where  the  Pqynting  vector  has  the  form 

£  =  X  £)  (18) 

and  the  volumetric  power  deposition  is  expressed  in  terms  of  the  total  electric  field  as 

p = ifleir  •  a-  (19) 

The  total  electric  field  £  in  the  geometric  optics  approximation  is  a  smn  over  individual 
geometric  optics  wave  fields 

£  =  E£>  (20) 

j 
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(21) 


and  the  total  magnetic  field  B.  is  likewise  decomposed, 

B.  = 

) 

For  each  geometric  optics  wave, 

Rj  =  -kj^Rj  (22) 

Ui 

and  the  total  current  density  Z  is 

i  =  (23) 

} 

In  order  to  evaluate  the  power  deposition  at  an  arbitrary  field  point,  the  conductivity 
tensor  must  also  be  interpolated  for  each  radiation  mode.  When  the  power  deposition 
is  evaluated  throughout  a  volume  in  the  ionosphere,  an  integral  form  of  the  energy 
conservation  equation  (17)  is  routinely  computed  and  chedced.  The  power  dissipated  in 
each  part  of  the  ionosphere  the  HF  radiation  is  compared  with  the  net  HF  power 
entering  that  portion  as  determined  by  integrating  the  Poynting  vector  over  appropriate 
surfaces. 

The  power  dissipation  in  (17)  is  the  mechanism  of  HF  ohmic  heating  in  the  ionosphere, 
and  is  used  as  a  source  term  in  the  energy  transport  analysis  discussed  in  Section  6. 

4-3.  Model  for  the  Background  Ionosphere 

As  already  mentioned,  the  geometric  optics  code  uses  a  modular  conductivity  tensor, 
that  is,  the  dielectric  properties  of  the  medium  are  incorporated  in  the  code  entirely 
through  an  isolated  module.  Of  the  various  ionosphere  modek  currently  in  usage(^) ,  we  have 
chosen  the  International  Reference  lonosphore  (IRl)  model^^^K  It  provides  electron  and 
ion  densities  and  temperatures  as  well  as  neutral  temperature  as  functions  of  geographic 
or  geomagnetic  location  and  altitude.  Parameters  are  time  of  day,  time  of  year,  and 
srinspot  number.  The  geomagnetic  field  is  also  evaluated  by  these  routines.  The  IRI  is  an 
empirical  model  that  gives  the  most  fiexibility  and  generality  for  applications  of  interest. 
The  nU  version  used  was  obtained  through  the  courtesy  of  D.  Bilitza  at  the  National  Space 
Science  Data  Center  at  NASA  Goddard  Space  Flight  Center.  The  neutral  density  model 
was  obtained  from  the  Mean  CIRA  Reference  Atmosphere  model^^^^  and  provides  densities 
of  molecular  nitrogen  and  oxygen,  atomic  oxygen,  argon,  helium,  and  hydrogen  at  various 
altitudes.  The  electron  collision  frequency  has  contributions  from  the  Coxilomb  collisions 
with  ions  as  well  as  short  range  colliaons  with  the  neutrals.  The  former  is  important 
in  the  F  region,  while  the  latter  dominates  the  D  and  E  regions  as  well  as  the  lower  F 
region  of  the  ionorohere.  The  magnetic  field  is  calculated  from  the  1968  POGO  spherical 
harmonic  modeli^^l.  Together  these  routines  provide  a  basis  for  the  determination  of  the 
electron  collision  frequency  and  dielectric  properties  of  the  ionosphere.  A  bicubic  spline 
interpolation  scheme  is  used  to  smooth  out  ixr^pilarities  in  the  IRI  electron  and  ion  density 
profiles.  This  procedure  ensures  continuous  fost  and  second  derivatives  for  the  electron 
density  and  temperature  profiles,  which  are  required  by  the  ray  tracing  code.  Figure  4-2 
shows  densities,  temperatures  and  cdlision  freciuency  profiles  used  in  our  computations 
and  derived  frxnn  the  IRI  routines. 
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Modular  design  of  the  software  allows  easy  substitution  of  models  for  properties  of  the 
ionosphere,  and  electron  density  in  the  ionosphere  can  also  be  modeled  with  an  analytical 
Chapman  profile.  Chapman  profiles  are  tailored  with  either  a  single  or  double  hump 
to  model  specific  nightime  (F  region)  or  daytime  (E  and  F  regions)  cases  of  interest. 
This  simpler  model  of  the  ionosphere  is  computationally  much  less  onerous  than  the  IRI 
model,  and  provides  a  convenient  means  to  quickly  predict  results  of  an  oblique  HF  heating 
experiment.  Parameters  in  the  Chapman  profile  may  be  determined  using  knowledge  of 
peak  density  and  its  altitude  and  scale  length,  usually  measured  during  experiments  with 
a  digisonde  or  a  backscatter  radar. 

4-4.  Other  Code  Improvements 

A  variety  of  minor  but  important  improvements  were  made  during  the  course  of  work  on 
other  tasks.  These  include  elimination  of  the  computation  of  \mimportant  dissipation  rate 
derivatives  for  the  ray  tracing  and  correction  of  an  inappropriate  form  of  the  vacuum  wave 
equation  tensor.  Spurious  trajectory  variation  terms  were  eliminated  from  the  plasma 
constitutive  relation,  and  the  definition  of  the  computed  phase  function  was  changed 
from  an  absolute  to  a  relative  one  for  numerical  reasons.  Two  coding  errors  in  the 
conductivity  routine  and  one  in  the  caustic  phase  jump  implementation  were  corrected, 
and  the  procedure  for  handling  change  of  degeneracy  of  the  background  was  modified  for 
convenience  in  subsequently  constructing  the  total  wave  field.  Source  code  documentation 
was  revised  and  improved. 

In  the  vicinity  of  caustic  surfaces,  where  radiation  modes  reverse  their  upward 
propagation  and  begin  to  return  to  earth,  the  field  amplitudes  are  locally  enhanced 
and  the  ionospheric  plasma  density  has  its  largest  values  along  the  propagation  path. 
Therefore,  the  strongest  radiation  ohmic  heating  effects  occur  in  these  regions  and  it 
is  partictilarly  important  for  self-consistent  radiation  propagation  calculations  to  have 
an  accurate  description  of  the  field  strengths  near  caustic  surfaces.  Unfortimately,  the 
geometric  optics  approximation  breaks  down  near  a  caustic  surface,  and  the  geometric 
optics  field  amplitude  actually  diverges  at  the  caustic  surface.  Presently,  the  geometric 
optics  code  computes  the  field  amplitude  using  a  technique  which  completely  avoids 
numerical  difi^culties  at  this  singularity  and  implements  a  jiunp  condition,  derived  from 
a  special  boimdary  layer  analysis,  for  the  geometric  optics  field  at  caustic  surfaces.  This 
allows  the  ray  tracing  calculation  to  proceed  accurately  in  the  region  beyond  the  occurrence 
of  the  caustic  singularity  but  does  not  give  accurate  values  of  the  field  amplitude  near  the 
caustic  surface.  This  defect  must  be  remedied  in  the  futiue. 

Other  improvements  to  be  made  in  the  future  include  the  following.  The  field 
interpolation  procedure  described  in  Section  5  must  be  improved  near  the  caustic,  and 
a  method  for  predicting  the  fields  in  the  shadow  region  near  the  caustic  must  be  devised. 
Also,  the  interpolation  algorithm  must  be  improved  in  order  to  handle  the  special  situations 
at  transitions  between  isotropic  (degenerate)  and  anisotropic  (nondegenerate)  media. 
In  order  to  accurately  predict  sp\irious  Doppler  shifts  due  to  time-dependent  plasma 
nonuniformities,  the  location  of  the  stirfaces  which  separate  degenerate  and  nondegenerate 
regions  must  be  precisely  determined.  A  simple  Newton  iteration  can  be  applied  when 
crossing  such  a  surface  to  solve  this  problem.  Re^Hcw&  uf  the  treatment  of  refiection  of 
radiation  from  the  earth  surface  and  of  the  progr.’tm  )<><;ic  associated  with  transitions 
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between  degenerate  and  nondegenerate  regions  are  planned.  If  calculations  of  nearly 
vertical  propagation  are  required,  an  improvement  in  the  ray  tracing  algorithm  is  called 
for;  however,  we  do  not  anticipate  such  reqtiirements  immediately  since  nearly  vertical 
propagation  is  of  little  interest  in  the  context  of  HF  radar  apart  from  the  modehng  of  HF 
ionospheric  sounding  by  peripheral  systems  in  an  HF  radar  installation.  A  special  case, 
which  has  never  yet  occurred  in  practice  but  may  occur  in  HF  skywave  calciilations,  is 
the  simultaneoxis  occtirrence  of  mode  separation  due  to  the  onset  of  plasma  anisotropy 
and  caustic  reflection.  The  ray  tracing  code  cannot  handle  this  event  and  it  may  become 
necessary  to  improve  the  code  accordingly. 

Figure  4-1.  HF  Antenna  Model 
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Figure  4*2.  Ionospheric  Parameters  Calculated  at  45° 
North  and  70°  West  for  the  Month  of  January  at 
8:00  AM  Local  Time  and  Sunspot  No.  100 
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5.  MODELING  OF  OBLIQUE  IONOSPHERIC  PROPAGATION 

A  model  has  been  developed  and  implemented  for  the  dielectric  tensor  of  the  ionosphere 
at  mid-latitude  and  high-latitude  using  standard  prescriptions  for  electron  density  and 
temperature  profiles. 

The  improved  wave  propagation  code  Has  been  exercised  to  scan  frequency, 
polarization,  latitude,  and  tilt  parameters,  which  affect  the  propagation  of  HF  beams 
in  the  ionosphere. 

In  this  section  we  show  some  illustrative  calculations  simulating  an  HF  r2ular  beam 
propagating  in  the  ionosphere.  All  examples  use  the  ionospheric  conditions  shown  in  Figure 
4-2.  Rays  are  transmitted  from  a  location  45°  North  and  70°  West  during  the  month  of 
January  at  8:00  AM  local  transmitter  time.  The  sunspot  number  is  100.  Figure  5-1  shows 
an  angular  scan  of  ten  rays  transmitted  at  10  MHz  with  horizontal  transmitter  polarization. 
Since  the  reflection  layer  altitude  is  a  function  of  frequency  and  angle  of  incidence,  some 
rays  reflect  from  the  ionosphere  and  others  pass  through  it. 

All  rays  bifurcate  into  distinct  ordinary  and  extraordinary  mode  rays  because  of  the 
anisotropy  of  the  magnetized  ionospheric  plasma.  Figure  5-2  shows  multiple  boimces  of  a 
radar-like  beam. 

Notice  that  while  rays  incident  at  steeper  angles  are  about  to  complete  four  hops, 
grazing  incidence  rays  are  barely  completing  two  hops.  A  full  layout  of  the  different 
projections  of  a  single  hop  beam  of  rays  is  shown  in  Figure  5-3. 

Details  of  the  caustic  surface  geometry,  splitting  of  the  beam  into  ordinary  and 
extraordinary  modes,  and  shifting  of  the  beam  pattern  toward  the  direction  of  the  earth’s 
magnetic  field  are  clearly  apparent.  Figure  5-4  shows  contours  of  constant  ohmic  power 
deposition  along  the  equatorial  x-z  plane. 

Here  the  contributions  from  different  waves  are  incoherently  superposed  since  the 
interference  terms  are  rapidly  oscillating  in  space.  Large  absorption  of  the  HF  energy 
occurs  at  low  altitude  in  the  D  and  E  regions  due  to  large  collision  frequencies.  Caustic 
swelling  at  higher  altitudes  is  also  a  cause  of  enhanced  absorption  of  the  wave  energy. 
These  regions  are  therefore  likely  locations  for  thermal  modification  of  the  ionosphere. 
The  code’s  capability  to  conserve  energy  is  carried  out  every  time  power  absorption  is 
calciilated,  in  the  manner  disciissed  in  Section  4. 

A  simple  way  to  test  how  density  irregularities  affect  the  radar’s  HF  beam  is  to  include 
a  density  blob  or  bubble  at  the  center  of  the  HF  beam.  Varying  the  amplitude  and  the  size 
of  the  irregrilarity  gives  a  parameter  regime  over  which  various  phenomena  are  likely  to 
occur.  Estimates  of  the  magnitude  of  effects  like  beam  self-focusing,  striations,  and  beam 
scattering  can  be  readily  obtained  in  this  non  self-consistent  manner. 

For  example,  in  attempting  to  find  out  how  large  a  density  blob  must  be  before  the 
path  of  an  initial  beam  of  rays  is  significantly  affected,  a  three  dimensional  gaussian  ball 
of  density  was  added.  Its  analytical  aq>ression  is  given  in  Section  8,  equation  (15),  and 
the  results  are  shown  in  Figures  8-9  to  8-11. 

Another  feattire  of  the  linear  analysis  of  oblique  propagation  is  the  modeling  of  ducted 
waves.  With  proper  laimch  angle  and  ionospheric  conditions,  rays  laimched  from  the 
grotmd  can  be  trapped  or  ducted  in  the  ionosphere  for  very  long  propagation  paths. 
Figure  5-5  shows  a  20  MHz  HF  1>cam  launched  from  the  east  coast  of  the  USA  (45°N, 
70°W)  at  an  angle  of  5.5°  above  horizontal  (other  parameters  are  same  as  before).  The 
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code  runs  show  that  the  rays  remain  ducted  for  over  6000  km  due  east  of  Moscow,  Maine 
and  could  potentially  remain  in  this  mode  for  longer  ranges,  i.e.,  without  hitting  the 
ground  or  escaping  from  the  upper  ionosphere.  Ducting  calculations  can  be  repeated  for 
other  operating  scenarios  amd  can  form  the  basis  of  a  separate  study  of  such  topics  as  the 
amoimt  of  HF  energy  ducted,  the  attenuation  suffered  by  a  ducted  beam,  induced  ducting 
by  by  powerful  HF  beams,  induced  duct  steering,  and  the  potential  for  reliable  detection 
of  reentry  vehicles  and  other  missiles. 

Figure  5-1.  Angular  Scan  of  Eastward  Propagating,  Horizontally 
Polarized  Rays  at  10  MHz;  Rays  Reflected  from  the  Ionosphere 
are  Reflected  Again  from  a  Curved  Earth  Surface  Assumed  to 

Be  a  Perfect  Conductor 
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Figure  5-4.  Power  Deposition 


p  -  <J .  E>  -  ^  J*-  E ,  with  E  -  ZEj  and  J  -  Z  J]  •  Ej 


Figure  5-5.  Propagation  of  20  MHz  Ducted  Beam  Launched 


6.  OVERVIEW  OF  THE  IONOSPHERIC  TRANSPORT  DESCRIPTION 

In  spite  of  increasing  interest  in  ionospheric  modification  physics  as  well  as  improved 
analytical  and  experimental  diagnostic  tools,  there  has  been  a  limited  number  of  attempts 
to  imderstand  and  simulate  bulk  processes  affecting  the  backgroimd  plasma  during  HF 
heating.  Plasma  transport  in  the  ionosphere  caused  by  enahnced  absoprtion  of  powerful  HF 
waves  is  one  such  important  effect.  Through  the  years,  it  has  been  modeled  an2dytically  and 
numerically  for  various  special  cases.  For  example,  mid-latitude  tranport  models  have  been 
developed  for  the  heating  experiments  conducted  at  Arecibo^^^’^^’^^l  and  Plattevillet^®’^®!. 
Transport  in  the  polar  ionosphere  has  also  been  modeled^^^,  as  has  the  case  of  a  colinear 
field-aligned  geometry,  for  which  the  effects  of  microwave  beam  heating  from  power  station 
satellitesf^^^  have  been  simulated  numerically.  Nmnerical  models  of  the  transport  equations 
have  been  derived  with  assumptions  which  either  limit  their  usage  to  a  specific  application 
or  which  are  simply  incomplete. 

In  the  current  work,  we  construct  a  transport  code  which  can  operate  equally  in 
different  regions  and  conditions  of  the  ionosphere.  The  current  transport  code  is  designed 
to  predict  temperature  and  density  variations  of  the  plasma  caused  by  HF  heating  in  the 
D,  E,  and  F  regions  of  the  ionosphere,  i.e.,  from  65  km  altitude  to  the  upper  ionosphere. 
It  operates  at  low,  middle,  and  high  latitudes.  It  includes  all  the  necessary  basic  atomic, 
molecular  and  chemical  ingredients  that  permit  a  valid  analysis  of  how  the  background 
reacts  to  ohmic  absoprtion  of  powerful  HF  waves.  The  transport  code  used  in  this  work 
was  developed  initially  to  model  high  latitude  polar  plasma  heating  in  the  F  region^^^. 
The  original  code  modeled  a  single  ionic  species  plasma  and  assumed  an  imperturbed 
equilibrium  backgroimd.  Figure  6-1  shows  an  example  of  results  from  the  polar  heating 
code  where  both  electron  temperature  and  density  are  changed  by  90  seconds  of  vertical 
heating  at  5  MHz. 

The  polar  transport  code  was  later  modified  to  include  the  effect  of  recombination 
on  F  region  heating^^^l.  It  is  this  later  version  of  the  code  that  has  been  upgraded 
to  include  the  lower  ionosphere  with  its  molecular  ions,  and  the  ability  to  operate  at 
arbitrary  geomagnetic  longitudes  and  latitudes.  Temperature  dependent  recombination 
rates  is  another  feature  added  to  the  current  code.  Section  7  discusses  the  details  of  the 
current  transport  model  as  well  as  its  numerical  implementation. 
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Figure  6>1.  Deformation  of  the  Plasma  Temperature  and 
Density  Under  90  Seconds  of  HF  Heating  Followed  by  60  Seconds 
of  Cooling.  HF  Signal  is  5  MHz,  EBP  =100  MW  (firom  Reference  6). 


(  I  SNi.  X4I00I 
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7.  UPGRADE  OF  THE  IONOSPHERIC  TRANSPORT  COMPUTER  CODE 

This  section  presents  the  transport  analysis  and  code  developments  which  have  been 
performed.  Previous  limitations  on  the  transport  calculation  boimdaries  have  been 
removed  so  that  transport  processes  in  the  lower  ionosphere,  below  150  km,  may  be 
modeled.  Phenomenology  of  the  E  region  has  been  incorporated.  The  molecular  ions 
NO'*'  and  Oj,  dominant  in  the  E  region,  have  been  included  in  the  model  along  with  O'*", 
which  dominates  the  F  region.  The  modeled  ionospheric  plasma  is  thus  composed  of  three 
ionic  species,  NO'*',  Oj,  and  O'*',  hereafter  denoted  as  species  1,  2,  and  3  respectively.  Also 
included  in  the  upgrade  are  the  effects  of  ionization,  recombination,  and  elastic  and  inelastic 
atomic  coUision  processes.  Although  densities  and  cross  sections  are  highly  tmeertain  in  this 
coUisional  region  of  the  ionosphere,  there  exists  sufficient  data  in  the  literature  to  permit 
qualitative  modeling  of  the  major  processes  occxiring  in  the  E  region.  The  transport  code 
has  been  generalized  to  accoimt  for  inclinations  of  the  geomagnetic  field  lines  at  arbitrary 
latitudinal  locations. 

Finally,  software  improvements  have  been  made  which  permit  more  flexible  interaction 
with  the  wave  propagation  code  as  well  as  with  other  routines  which  model  background 
ionospheric  parameters. 

7-1.  Momentum  Equations 

The  momenttun  equations  for  the  four  species  (3  ionic  +  1  electronic)  are  derived  from 
fitiid  theory,  by  making  the  following  assumptions: 

ions  are  singly  ionized  (Z  =  1);  (1) 

ni  +  n2  +  ^3  =  ^<  »  (quasi  —  neutrality);  (2) 

+  n23i2  +  ”3Ji3  =  »  (ambipolar  diffusion);  (3) 

where  is  the  fluid  velocity  of  species  a.  When  inertia  is  neglected,  the  momentum 
equations  take  the  form 

Vpo  =  naqaK  +  -na^alio  ^  R-  maTlaUa^an  “  rnaUa  (tio  “  U.0)^a0  ~  ,  (4) 

where  a  =  1, 2,3,4  refers  to  the  ionic  species  1,2,3,  and  the  electrons,  respectively.  The 
coUisional  rate  Van  represents  collisions  of  species  a  with  the  neutral  constituents,  which 
are  atomic  hydrogen  (fr"),  atomic  helitun  {He),  atomic  oxygen  (O),  molecular  nitrogen 
(N2),  and  molecular  oxygen  (02)-  The  coUisional  rate  1/^0  represents  electron-ion  and 
ion-ion  coUisions.  The  quantity  go  is  the  charge  of  species  a;  ttiq  is  the  mass  of  species  a; 
c  is  the  speed  of  light,  and  the  pressure  of  species  a  is 

Pa  =  UqTq  .  (5) 

Acceleration  due  to  gravity  is  denoted  by  £. 

When  equation  (4)  is  summed  over  the  species  index  a,  then 

4 

V(pe+pi+P2+P3]  =  -menet'enlie-minit/ittJll-m2n2U2,i3i2-^3^3^nli3-(J2  ’^o"a)£  ■ 

a=l 

(6) 
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The  terms  proportional  to  the  electric  field  £.  vanish  by  invoking  quasi-neutrality,  and  those 
proportional  to  vanish  by  invoking  ambipolar  diffusion.  The  terms  proportional  to 

the  electron-ion  and  ion-ion  coUisional  rates  cancel  exactly. 

Since  this  model  supports  ambipolar  diffusion,  it  is  necessary  to  determine  the 
fractional  fliix  for  each  ionic  species  given  a  unit  fiux  of  electrons.  Such  a  quantity  depends 
upon  the  mass,  density,  and  coUisional  rate  of  each  species.  Thus,  it  is  assumed  that 


^aiia 

where  a  s  1,  2,  3  (ionic  species  index), 


c  MlUiff 
—  Oa 

^a^an 


c  — 

Oq  =  , 

ne 


(7) 


(8) 


and  Mjvjff  is  a  mass-weighted  ionic  coUisional  term,  which  is  determined  by  summing  (7) 
over  a  and  using  (3), 


^I^IN  = 


6i(m2U2n)im3t^n)  + 
Substitution  of  (7)  into  (6)  yields 

-1 


rieJU  ~ 


•  { V(p«  4-  PI  +  P2  +  P3]  +  E  mana)£}  . 


(9) 


(10) 


rrieUen  +  Mjujj\r  a 

When  the  projection  of  the  electron  fiux  onto  a  magnetic  field  line  is  taken,  it  is  foimd  that 


neWe  s=  -U(— (ne(T«  + 1;)]  -f-  '^anag  •  R/B]  , 


(11) 


where  Wa  =  Ua  '  R/B,  s  =  r  •  R/B,  R  =  Bxix  +  Byty  -I-  Bg^g,  B  =  [B|  -I-  Sj  -|-  is 

the  earth’s  magnetic  field,  and  the  diffusion  coefficient  D  is 


D  = 


(12) 


and  where  it  has  been  assiuned  that  the  ionic  species  are  at  the  same  temperature,  namely, 

Ti  =  Ti  =  T2  =  Ts  .  (13) 

The  ion  (7^)  and  neutral  (Tn)  temperatures  are  treated  as  constant  in  time,  but  are 
aUowed  to  vary  in  space,  a  feature  which  was  not  included  in  past  versions  of  the  transport 
code(^’^^).  Equations  (10),  (11),  and  (12)  assume  the  ambipolar  plasma  form  with  a  single 
ion  species  having  an  effective  mass  of  Mj  and  an  effective  ion-neutral  coUision  frequency 
oi  vjff.  A  stimmary  of  the  coUisional  rates  included  in  this  model  is  given  in  Appendix  A. 

7-2.  Continuity  Equations 

The  continuity  equation  along  a  magnetic  field  line  for  each  species  is  of  the  form 


dn  d 

'W  ^  ~  )®o(«,<)  +  Pai^^t)  +  Soa{s)  , 


(14) 
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where  Soa  is  the  steady  state  source  term  which  brings  about  the  equilibrium  density 
profile,  Pa  is  the  production  rate,  and  $a  is  the  recombination  rate  for  species  a. 

According  to  Schunk  and  Walker^^®^ ,  the  important  chemical  reactions  from  100  —  600 
km  in  altitude  which  involve  O'*",  and  NO"^  are 


0++JV2 iVO++JV  , 


k\  =  12- 


Tn  -  300  .  (Tn  -  aOOKTn  -  600)  .,,.13 
50  3.5  X 10^ 


cm^/s  (15) 


0+  +  O2  Of  +  O  ,  ki  =  1.1  X  10-®/!?  '^  cmVs  (16) 

0^  +  N2-^  NO"^  +  O  ,  *4  =  5.0  X  10"^*  cm*/8  (17) 

JVO+  +  c"  AT  +  O  ,  k'lo  =  4.2  x  10"'^(300/Te)  ®®  cm^s  (18) 

Of  +  c"  O  +  O  ,  k\i  =  1.6  X  10"‘^(300/re)  “  cm^/s  .  (19) 


The  dissociative  reaction  rates  kjQ  and  kji  are  those  cited  by  Schunk  and  Nagy^^®^  The 
reaction  rates  kj  and  k^  are  those  cited  by  Park  and  Banks^^®].  All  other  reaction  rates 

are  from  Schunk  and  Walker^^®!. 

The  continuity  equations  then  take  the  form 


— ni  +  — (nitui)  =  Sol  +  ’»(^2)[*W3  +  *'4»»2]  “  (20) 

Q  Q 

^rx2  +  ^(n2tr2)  =  <^o2  +  ^p^{02)nz  -  [*;h»»e  +  k4n(iV2)]n2  (21) 

=  •^03“  I*:i«(-^2)  +  ^»»(^^2)]^3  (22) 

Q  Q 

^ne  +  ^(«etUe)  =  5oe  -  [kjoni  +  k'iin2]ne  .  (23) 


The  molecular  ions,  ni  and  n2  are  assumed  to  be  in  chemical  equilibrivun,  in  accord  with 
the  findings  of  Schunk  and  Walker^^®).  Hence,  the  lefthand  sides  of  equations  (20)  and  (21) 
are  set  to  zero,  and  these  equations  then  reduce  to 


Sol  +  {M»»3  +  *4”2]»’(<^2) 

_  So2  +  k2n(  02)113 

~  +  M(<^2) 

The  quasi-neutrality  equation  (2)  implies  that  the  sum  of  (24)  and  (25)  yields 


na  =  Tie  -  ni  -  ^2 


klone  +  k;n(JV2)  -h  k^n(02)^^^^^ 


(24) 

(25) 


(26) 


Originally,  Shoucri  et  al.(^)  did  not  include  recombination  in  the  continuity  equation, 
as  it  was  their  desire  to  model  short  time  {t  <  90.<)  effects  in  the  F  region.  They  found 
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that  the  density  changes  caused  by  the  pressure  imbalance  induced  by  HF  heating  never 
sattirated  over  these  short  time  periods.  Later,  Hansen^^^^  added  a  recombination  rate  to 
the  continuity  equation  of  the  form 

r„)  =  4n(02)  +  *'in(JV2)  s'*  .  (27) 

This  recombination  rate  modds  dissociative  recombination  of  NO"^  and  by  a  loss  of 
O'**,  an  effect  which  is  currently  present  in  the  continuity  equation  (22)  for  O'**  as  well  as 
in  the  equations  (20)  and  (21)  for  NO"^  and  as  production  rates. 

The  reaction  rates  for  NO"^  and  on  functions  of  the  electron  temperature  Te,  and 
as  indicated  in  equation  (23),  there  is  a  loss  of  electrons  at  a  rate  which  depends  upon  the 
electron  temperature  in  a  nonlinear  fashion,  namely,  fitiTt)-  Physically,  what  occurs  is 
that  the  ionosphere  is  heated,  and  the  electron  temperature  increases;  the  electrons  then 
do  not  recombine  as  easily  as  they  did  before  heating.  This  is  modeled  in  (18)  and  (19) 
by  the  reaction  rates  being  inversely  proportional  to  the  electron  temperature  raised  to  a 
power  which  has  been  determined  empirically^^®!. 


7-3.  Electron  Energy  Equation 

The  electron  energy  eqtiation  along  a  magnetic  field  line^^!  is  of  the  form 

+  =  +  +  ,  (28) 

where  k  is  Boltzmann’s  constant,  and  where  the  first  term  represents  the  time  rate  of  change 
of  the  electron  temperature  Te  due  to  heating;  the  second  term  represents  convection,  and 
the  third  pressure  fiux.  The  first  term  on  the  righthand  side  of  (28)  represents  thermal 
conductivity  along  a  magnetic  field  line,  where  Ke  is  the  parallel  coefficient  of  thermal 
conductivity!®!, 


K« 


7.75  X  10®2f 

1  +  3.22  X  <Qk> 


eV  cm“^  s”^  K  ^ 


(29) 


and  where  the  summation  in  the  denominator  is  over  all  the  products  of  neutral  gas  species 
concentrations  with  average  electron-neutral  momentum  transfer  cross  sections  <  Qjg  >. 

The  second  term  on  the  righthand  side  of  (28),  Q],f,  is  the  HF  heating  source, 
Qlif  =  RefEJj  •  J]  =  Re[l^  •  £  •  where  fiijf  is  the  complex  amplitude  of  the  HF 
electric  field  whose  full  space-time  dependence  is  given  by  [Ehf(t)e~“*^*  +  c.c.]/2.  The 
third  term  on  the  righthand  side  of  (28),  Qo,  accounts  for  the  steady  state  source  which 
brings  about  the  equilibriiun  temperature  profile  in  the  absence  of  the  HF  pump.  The 
last  term  on  the  righthand  side  of  (28),  represents  the  rate  of  electron  cooling  caused 
by  translational,  rotational,  and  vibrational  excitation  of  neutrals  as  well  as  cooling  by 
collisions  with  ions.  A  sxunmary  of  the  terms  included  in  Le  can  be  found  in  Appendix  B. 


7-4.  Recovery  of  Limiting  Cases 


In  the  limit  of  ni  s  n2  ~  0  the  transp<»t  eqtiations  used  by  Hansen!^®!  are  recovered. 


d 


Q 

|jj(n3U>3)  =  5<,3  -  (fcin(7\r2)  -|-  i:2n(02)ln3 


(30) 


30 


and 


rieWe  =  -D(^[ne(Te  +  T;)]  +  ^a^ag  ’  R/B]  ,  (31) 

where  >13  =  n*  (qviasi-neutrality),  and  D  =  electron  energy  equation 

used  by  Hansen^^^I  is  the  same  as  (28),  except  that  some  updates  have  been  made  to  the 
Le  term  for  the  current  model. 

If  the  gravitational  term  and  the  recombination  term,  which  are  terms  linear  in  the 
electron  density  or  its  first  derivative,  are  incorporated  into  the  equilibrium  source  term 
in  (30),  an  assiunption  valid  when  heating  on  short  time  scales,  then  the  transport  model 
of  Shoucri  et  al.^®!  is  recovered, 

Q  a 

+  ^(n«u>e)  =  Soe  (32) 

= -i)l^(ne(Te  +  ri)l  .  (33) 

The  electron  energy  equation  used  by  Shoucri  et  al.^®^  is  identical  to  that  of  Hansen^^^^ 
The  ciurent  transport  model  has  the  ability  to  recover  these  limits.  It  can  be  run  with 
the  recombination  rate  updated  in  time  (because  Te  is  changing),  it  can  be  run  with  the 
recombination  rate  constant  in  time  (in  analogy  to  Hansent^^]),  or  it  can  be  nm  without 
the  effect  of  recombination  (Shoucri  et  al.(®l) 

7-5.  Numerical  Solution  for  Transport 

Using  a  predictor-corrector  form  of  the  time-averaged  Crank-Nicholson^^^^  tri-diagonal 
method,  eqtiations  (11),  (12),  (23),  (24),  (25),  (26)  and  (28)  are  solved  when  the  HF  energy 
is  turned  on.  The  temporal  and  spatial  modifications  of  density  and  temperatme  are  solved 
in  the  following  manner: 

1.  Solve  for  T®  using  (28)  with  Qhf  turned  on. 

2.  Update  all  coefficients:  Le,  Ke,  i/en* 

3.  Solve  for  the  electron  density,  Ue,  using  (23). 

4.  Update  all  coefficients:  Le,  Ke,  Uta, 

5.  Update  the  O'*"  density  (n3),  the  density  (*12),  and  the  NO'^  density,  using  (24), 
(25),  and  (26). 

6.  Update  all  coefficients:  Le,  Ke,  i^eni 

7.  Repeat  at  each  time  step. 

The  botmdary  conditions  utilized  are  those  established  by  Hansen^^^K  Any  linear 
combination  of  Te  and  dTe/ds,  or  ne  and  driefds,  may  be  used  as  boimdary  conditions  at 
the  lower  and  upper  endpoints  for  equations  (23)  and  (28).  For  the  current  nms,  at  the 
lower  boundary. 


re(-  —  ~  T«{z  —  ZiQy,gf.,t  —  0) 
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(34) 

(35) 


whereas,  at  the  upper  botmdary, 

dT^fds{z  —  ^upperi^)  —  &Tflds{^z  =  z^pper^'t  =  0)  (36) 

dne/ds{z  =  Zupperti)  —  dnefds^z  =  =  0)  (37) 

In  the  D-E  layer,  it  has  been  determined  that  the  temperature  and  density 
modifications  are  essentially  local  in  nature  because  the  thermal  conductivity,  compression, 
convection,  and  diffusion  coefficients  nearly  vanish.  Thus,  in  this  region  of  space,  the 
conductivity,  compression,  and  convection  terms  are  set  equal  to  zero  in  (refte  eqn),  and 
the  local  balance  between  heating  and  cooling  determines  the  temperature  modification, 
which  is  isotropic  in  nature,  because  Uce  In  the  F  layer,  however,  where  these  terms 
are  not  negligible,  the  complete  equation  (28)  is  solved,  and  the  temperature  modifications 
are  field-aligned  because  cjce  »  Vf 

7-6.  Establishment  of  the  Equilibrium 

Before  any  heating  can  be  done,  it  is  necessary  to  first  set  up  the  equilibrium  profiles. 
For  Te,  Tj,  Tn,  n^,  and  the  ionic  densities,  this  is  carried  out  by  invoking  IRI  (International 
Beference  Ionosphere)  at  the  latmch  point  to  obtain  a  data  set,  and  by  then  making  a 
cubic  spline  fit  to  the  data  as  a  function  of  altitude  only.  This  guarantees  continuous 
first  and  second  derivatives,  necessary  for  both  the  ray  tracing  and  transport  codes.  While 
variations  of  temperature  and  density  with  respect  to  latitude  and  longitude  have  not  been 
implemented,  the  capacity  for  including  such  variations  has  been  developed. 

The  neutral  densities  are  set  up  using  CIRA  data,  and  there  is  not  a  spline  fit 
incorporated  into  this  package  for  these  profiles.  Since  there  is  scarce  CIRA  data  for 
atomic  hydrogen  density,  it  was  modeled  using  a  Boltzmann  law  for  the  pressure  as  a 
function  of  altitude. 

Two  typical  cases  have  been  run  using  the  transport  code.  The  first  is  due  east  HF 
propagation  &om  Moscow,  Maine,  with  stinspot  num^r  150  (a  value  typical  for  1991  and 
1992),  in  January  at  1200  hours  local  time,  denoted  as  case  0112.  Depicted  in  Figtnres  7-1, 
7-2,  and  7-3  are  the  neutral  density  profiles;  the  equilibrium  profiles  of  Te,  7J;,  and  Tn;  and 
the  equilibrium  profiles  of  ne,  nj,  n2,  and  for  case  0112. 

The  second  case  is  due  east  propagation  from  Moscow,  Maine,  with  sunspot  ntimber 
150  in  July  at  0400  hours  local  time,  denoted  as  case  0704.  Figures  7-4  and  7-5  exhibit  the 
equilibrium  temperature  and  density  profiles  for  this  case.  (The  neutral  density  profiles 
are  the  same  as  in  case  0112.) 

The  original  transport  codet®!  was  developed  for  the  polar  region  and  hence  Vn  H  £  |1 
for  such  geometry.  Later,  Hansenl^^l  extrapolated  the  code  to  describe  mid-latitude  effects 
such  as  those  observed  at  Aredbo  Observatory,  and  th\is  the  magnetic  field  was  chosen 
to  have  a  dip  angle  with  a  typical  value  of  50®.  In  the  current  version  of  the  code,  the 
transport  calcxilation  is  run  on  a  magnetic  field  line  whose  orientation  is  specified  as  an 
input.  Given  the  beginning  point  of  a  field  line  and  the  grid  spacing,  the  field  line  along 
which  modifications  are  to  be  calctilated  is  established.  Associated  with  each  grid  point 
on  the  field  line  are  two  reference  frames:  the  latmch  point  reference  frame  (l.p.),  where 
(x,y,z)  =  (0,0,0)  is  at  the  launch  point;  and  the  center-of-earth  reference  frame  (c.o.e.), 
which  provides  the  latitude,  longitude,  and  altitude  of  each  grid  point. 
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Once  an  altitude  is  associated  with  each  grid  point  on  a  field  line,  the  eqtiilibrium 
densities  and  temperatiires  of  that  grid  point  can  be  specified.  Latitudinal  and  longitudinal 
variations  in  the  density  and  temperature  profiles  can  be  included  by  calling  IRI  data  for 
each  grid  point. 

It  is  necessary  that  the  first  and  second  derivatives  of  density  and  temperature  be 
continuous  for  the  transport  calculations,  because  the  modifications  of  temperature  and 
density  depend  upon  their  derivatives,  which  are  incorporated  into  the  diffusion  term,  the 
thermal  conductivity  term,  the  convective  term,  and  the  pressxure  fiux  term.  The  cubic 
spline  fit  to  the  IRI  data  addresses  this  requirement. 

When  the  equilibrium  profiles  are  established,  the  sources  which  sustain  the  profiles 
are  determined.  Equations  (23)  through  (26)  are  solved  in  steady  state  to  ascertain  the 
density  sources,  Soa-  The  temperature  source  Qo  is  foimd  from  (28)  with  d/dt  =  0  and 

Qhf  =  0- 


7-7.  Versions  of  the  Transport  Code 

Three  different  versions  of  the  transport  code  currently  exist.  The  first  is  “eos-oth”, 
which  runs  the  transport  on  many  field  lines  in  many  different  magnetic  meridian  planes. 
The  second  is  “eos_oth_lplane”,  which  runs  the  trzmsport  on  many  field  lines  in  just 
one  magnetic  meridian  plane.  To  append  all  of  the  output  files  together  from  running 
“eos_oth_lplane”  on  many  different  planes,  there  exists  a  short  program  “output”.  Finally, 
there  is  a  version  of  the  code  called  “eos_diag”,  which  allows  the  operator  to  choose  a 
specific  field  line  in  a  specific  magnetic  meridian  plane  upon  which  he/she  desires  to  see 
diagnostic  output  of  all  equilibrium  quantitites,  as  well  as  the  quantities  which  vary  in 
time  because  of  heating.  The  program  “eosplot.oth”  plots  the  diagnostic  output  ^m 
“eos-diag”. 
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7-1.  Neutral  Density  Profiles 


Figure 
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Equlllbrlua  T*ap«ratur«M, 


Figure  7-2.  Ek|uilibrium  Temperature  Profiles  for  Case  0112 
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Figure  7-4.  Ekiuilibrium  Temperature  Profiles  for  Case  0704 
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Figure  7-5.  Equilibrium  Charged  Particle  Density  Profiles  for  Case  0704 
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8.  COUPLING  OF  THE  WAVE  PROPAGATION  AND  TRANSPORT 
CODES 


The  block  diagram  in  Figure  3-3  shows  the  coupling  of  the  transport  code  to  the  wave 
propagation  code.  Here  the  \ise  of  the  output  of  the  wave  propagation  code  is  used  as 
input  to  the  transport  code. 

An  algorithm  for  projecting  the  propagating  geometric  optics  rays  associated  with 
an  HF  beam  onto  the  field-aligned  transport  calculation  numerical  grid  has  been  devised 
and  implemented.  Initially,  the  rays  are  traced  through  the  ionosphere  to  determine  the 
imperturbed  trajectories.  Then  the  region  of  the  ionosphere  in  which  the  modification  is 
to  be  calculated  is  earmarked,  and  is  contained  within  a  voliune  spanned  by  the  points  rx , 
I2»  I3»  I4*  The  point  n  is  typically  chosen  so  that  it  is  below  the  bottom-most  ray  in 

altitude,  near  the  entry  point  to  the  ionosphere,  and  outside  of  the  beam  in  the  y-z  plane. 
The  point  12  is  then  chosen  so  that  the  vector  12  —  Ii  lies  along  a  magnetic  field  line,  with 
12  high  enough  in  altitude  that  all  of  the  HF  beam  is  contained  on  the  field  line. 

The  points  ^  and  14  cannot  be  chosen  until  the  magnetic  field  coordinate  system  is 
established.  A  magnetic  field  line  is  constrained  to  lie  along  the  direction  of  the  unit  vector 
63,  where 

Bx  By  Bz 

in  the  laimch  point  reference  frame,  and 


(1) 


63  =  (0, 0, 1)  , 

in  the  magnetic  field  coordinate  frame,  with 


(2) 


(3) 


The  unit  vector  52,  which  denotes  the  direction  to  adjacent  magnetic  field  lines  in  the  same 
magnetic  meridian  plane,  is 


S2  = 


BxBz  ,  ByBz 

-  — s 


BB^^  ’  BB^ 

in  the  latmch  point  frame,  and 

62  =(0,1,0)  , 

in  the  magnetic  field  coordinate  frame,  with 


Ea. 

B  ^  ’ 


5x=[s2  +  b2]1/2 


(4) 

(5) 

(6) 


Finally,  the  unit  vector  gx  defines  the  direction  to  an  adjacent  magnetic  meridian  plane, 
and  is  given  by  the  expressions 


_  By  Bx 


(7) 


in  the  launch  point  reference  frame,  and 


gl=  (1,0,0)  . 
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(8) 


in  the  magnetic  field  coordinate  frame.  There  is  a  well-defined,  unitazy  transformation 
which  links  the  two  reference  frames,  namely. 


where 


(”). 

\^/  B  frami' 


=  BT^ 


By/Bj_  B^BzIBB^  B^/B' 
B  =  I  ByBzfBB^  By/B 

0  -BJB  Bz/B^ 


(9) 


(10) 


(11) 


and  £—1  =  so  that  £  •  ^  =  1. 

The  point  £3  is  constrained  to  lie  in  the  plane  determined  by  the  points  r^,  -f  As3e3, 
Mid  Cl  +  A82fi2  >  the  g2  —  S3  plane,  which  is  that  of  the  magnetic  meridian.  Typic^y, 
£3  is  at  the  same  altitude  as  £i  (on  the  bottomside  of  the  beam  in  the  x-y  plane),  but  on 
the  side  of  the  beam  opposite  to  £]^  in  the  y-z  plane. 

Finally,  the  point  £4  must  lie  in  the  §3  —  S3  plane,  so  that 


(14 -Cl)  •  fi2  =  0  ,  (12) 

i.e.,  (C4  ~  £l)  has  a  component  across  meridian  planes  as  well  as  along  magnetic  field  lines, 
but  not  a  component  in  the  direction  between  adjacent  field  lin«^  in  the  same  meridian 
plane.  Figures  8-1,  8-2,  and  8-3  show  the  xmperturbed  rays  in  the  x-z,  y-z,  and  x-y  planes, 
respectively,  as  well  as  the  locations  of  £i,  £3,  £3,  and  £4  for  case  0112.  Figures  8-4,  8-5, 
and  8-6  depict  the  unperturbed  rays  in  the  x-z,  y-z,  and  x-y  planes,  respectively,  with  the 
locations  of  £x,  £3,  £3,  and  £4  denoted,  for  case  0704. 

Once  the  points  £1,  £3,  £3,  and  £4  are  determined,  a  numerical  grid  is  established. 
A  numerical  interpolation  of  the  electric  field  £  is  performed  at  specific  points  along  a 
magnetic  field  line,  whose  bottommost  grid  point  is 


So 


II  + 


14 -II 
#  planes  —  1 


X  i 


+ 


13-11 

#  lines  —  1 


X  j 


(13) 


where  i  <  #  planes,  and  j  <  #  lines.  The  spline  fit  call  is  made  to  retrieve  the  election 
density  and  temperature  at  the  point  of  interpolation  so  that  the  quantity 


Ohf  =  Re<£*  •  j£>.  (14) 

at  each  grid  point  along  each  field  line  in  each  meridian  plane  is  determined.  Samples  of 
Olif  profiles  for  cases  0112  and  0704  are  exhibited  in  Figures  8-7  and  8-8. 

An  effort  has  been  made  to  assess  just  how  large  the  density  modifications  mr^  be  in 
order  for  the  rays  to  be  affected.  Thus,  a  Gaussian  ball  of  density  was  non-self-consistently 
introduced  into  the  ray  tracing  code,  so  that 

n(l)  =  no(l)  +  exp{-(zV«^  +  y W  +  *^/c^)]  • 
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(15) 


The  parameter  can  be  either  positive  or  negative.  It  was  found  that  for  case  0112,  when 
6n  ~  .2no  and  is  situated  at  an  altitude  of  100  km  along  the  center  of  the  beam,  the  rays 
are  nominally  disturbed.  The  ray  paths  are  displayed  in  Figures  8*9  through  8-11  when 
this  density  modification  is  included. 

This  is  tantamoimt  to  allowing  the  density  and  temperature  modifications  to  satmate 
in  the  lower  ionosphere,  which  occurs  within  two  minutes  of  heating.  Based  upon  this 
analysis,  it  was  decided  to  allow  the  density  and  temperature  modifications  to  evolve  for 
two  minutes  before  performing  the  feedback  calculation. 

8>1.  Daytime  Heating  -  Case  0112 

For  the  daytime  scenario  described  in  Section  7-6,  Figures  8-12  through  8-18  display 
the  modification  which  occurs.  The  data  shown  is  firom  field  line  #6  on  plane  10,  where 
the  pump  is  centered  at  an  altitude  of  roughly  100  km.  Since  the  daytime  density  is  so 
high  2  X  lO^cm"^  at  the  peak  of  the  profile),  the  rays  never  get  higher  than  18^  km  in 
altitude.  The  frequency  of  the  HF  beam  is  15  MHz  for  this  case. 

Figure  8-12  depicts  the  coefficient  of  thermal  conductivity,  Ke,  normalized  to  aL^kncoo) 
where 

Lq  =  ionospheric  scale  length  =  5  x  10^ cm  (16) 

k  =  Boltzmann^s  constant  (17) 

n«oo  =  no(s  =  Supper)  (18) 

_  2kTeoo/l^ 

"  = 

Teoo  =  To(s  =  Supper)  (20) 

P  =  a  characteristic  ionospheric  collision  frequency  =  4.444  x  10^s“^  .  (21) 

Essentially,  Ke  is  zero  in  the  vicinity  of  the  beam.  Corrections  introduced  by  the  convective 
and  pressure  flux  terms  are  even  smaller  than  that  of  the  thermal  conductivity  term^^^. 
Thus,  for  this  case,  the  heating  is  local  in  nature,  as  is  evidenced  in  Figures  8-13a  and 
8-136.  Figure  8-13a  depicts  Te  vs.  altitude,  plotted  at  6  second  intervals  for  the  duration 
of  heating  (2  minutes).  Figiue  8-136  is  a  plot  of  6Te(s)  =  [Te(s,t)/Te(s,t  =  0)  —  1]  vs. 
altitude  at  6  second  intervals.  The  peak  modification  normalized  to  the  background  is  on 
the  order  of  75%,  and  is  self-similar  to  the  pump  shape.  A  typical  run  was  also  done  which 
included  the  thermal  conductivity,  convective  term,  and  pressure  flux  term.  There  is  no 
significant  change  in  the  output  when  these  terms  are  excluded  below  140  km  in  altitude. 
They  are,  however,  necessary  and  important  above  140  km  in  altitude,  and  are  retained 
when  the  ptunp  center  falls  above  this  point. 

Figure  8-13c  exhibits  6Te(s)  near  the  peak  of  the  profile  vs.  time.  The  changes  in  the 
electron  temperature  sattirate  quickly  on  a  short  time  scale;  there  is  also  a  slight  decrease  in 
6Te  on  a  longer  time  scale  ('^  minutes),  which  is  in  response  to  the  density  modifications. 
The  density  increases  on  the  long  time  scale  in  the  vicinity  of  the  pump,  and  thus  there 
are  more  particles  to  heat  than  in  the  early  heating  stage.  As  a  consequence,  the  electron 
temp>erature  slightly  decreases  on  the  longer  time  scale. 

Figxires  8-14  and  8-15  are  plots  of  the  electron  cooling,  Le(s)  vs.  altitude,  and  the 
electron  collision  frequency,  */e(s)  =  Um  +  J'ei  vs.  altitude,  respectively,  at  6  second 
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intervals.  The  electron  cooling  shows  an  increase  centered  around  the  ptunp  which 
satiurates  on  the  temperature  time  scale.  The  collision  frequency  also  shows  an  increase 
centered  about  the  pump.  However,  as  the  input  energy  moves  higher  in  altitude  on  a 
field  line  and  enters  the  region  where  the  collisions  are  electron-ion  in  nature  rather  than 
electron-neutral,  the  change  in  the  collision  frequency  becomes  a  decrement  rather  than 

an  increment,  because  is  inversely  proportional  to  Te  ,  whereas  i/en  scales  as  T®'  . 

Figure  8- 16a  depicts  n«(s)  vs.  altitude  plotted  at  6  second  intervals  over  the  course 
of  2  minutes  of  heating.  The  density  modification  normalized  to  the  background  density, 
^ne(s)  =  [ne(s,  t)/ne(s,  t  =  0)  —  1]  vs.  altitude  at  6  second  intervals  is  displayed  in  Figure 
8-166.  Figure  8- 16c  shows  6ne(s)  near  the  center  of  the  pump  as  a  function  of  heating  time. 
The  density  modification  is  self-similar  to  the  pump,  is  an  increase  rather  than  a  decrease, 
and  saturates  on  a  time  scale  of  minutes  rather  than  seconds.  The  density  increase  may  be 
explained  by  examining  Figures  8-17  and  8-18.  Figure  8-17  is  a  plot  of  the  recombination 
rate,  i6e(s),  vs.  altitude  at  6  second  intervals.  It  is  readily  observable  that  fit  decreases 
aroimd  the  pump  as  heating  occurs,  because  T^  increases  [c.f.  equations  (17),  (18),  and  (22) 
of  Section  7].  Figure  8-18  exhibits  the  diffusion  coefficient,  D(8),  normalized  to  In  the 
vicinity  of  the  pump,  D  is  essentially  zero,  and  the  recombination  rate  fit  decreases  from 
what  it  was  in  equilibrium.  However,  the  source,  Se,  is  emitting  enough  photoelectrons 
to  sustain  the  equilibrium  profile.  Thercdbre,  the  decrease  in  the  recombination  rate,  fit, 
caused  by  heating,  causes  an  increase  in  the  electron  density. 

The  features  of  the  modification  summarized  in  Figures  8-12  through  8-18  hold  on 
each  field  line  in  each  meridian  plane  wherever  modification  occurs  for  case  0112.  Only  the 
amplitudes  of  the  modifications  are  altered  from  field  line  to  field  line,  and  from  meridian 
plane  to  meridian  plane. 

8-2.  Nighttime  Heating  -  Case  0704 

Figures  8-19  through  8-25  describe  the  density  and  temperattire  modifications  for  a 
nighttime  scenario.  The  data  shown  is  from  field  line  #9  on  plane  20,  where  the  pump 
is  roughly  centered  at  an  altitude  of  170  km.  The  quantities  n«oo  aiid  Teoo  are  now 
normalized  to  the  lowest  grid  point  on  the  field  line  rather  than  the  highest  grid  point, 

neoo  =  ne(s  =  Slower)  (22) 

Teoo  —  Te(8  =  Slower)  (23) 

The  electron  density  peaks  at  a  value  of  about  3  x  10^ cm”^.  The  firequency  of  the  HF 
beam  is  set  to  10  MHz.  Figure  8-19  depicts  the  coefficient  of  thermal  conductivity,  Ke, 
normalized  to  aLoneoc*  the  vicinity  of  an  altitude  of  170  km,  Ke  is  non-zero,  and  thus 
it  is  anticipated  that  conduction  of  heat  away  from  where  the  pump  is  located  will  occur. 
For  these  nightime  runs,  all  of  the  terms  in  equation  (28)  of  Section  7  were  retained  when 
the  pump  center  was  higher  that  140  km  in  altitude. 

Figures  8-20  exhibit  the  temperatiue  modifications  which  take  place  for  case  0704. 
Figiire  8- 20a  is  a  plot  of  Te  along  a  field  line  vs.  altitude,  plotted  at  6  second  intervals 
over  the  course  of  2  minutes  of  heating.  There  is  clearly  a  change  in  Te  at  locations  other 
than  where  the  pvimp  heats.  Figure  8-206,  which  displays  6Te(s)  vs.  altitude,  shows  this 
phenomenon  more  clearly.  There  is  a  definite  conduction  of  heat  to  hiirlifr  altitudes.  The 
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change  in  ^Te  at  lower  altitude  is  not  as  great  because  Ke  is  vanishingly  small  there.  Figure 
8-20c  is  a  plot  of  ^Te  near  the  peak  of  modification  vs.  time.  The  electron  temperature 
response  time  to  the  heating  is  on  the  order  of  a  few  seconds. 

Figure  8-21  depicts  the  electron  cooling  term,  Le(s),  vs.  altitude  at  6  second  intervals. 
There  is  an  increase  in  Le  of  roughly  25%  around  the  pump  center,  and  it  essentially 
occurs  within  the  first  6  seconds.  Figure*  S-22  displays  the  electron  collision  frequency 
Ue(s)  =  i/en  +  ^ei  altitude  at  6  second  intervals  throughout  the  course  of  heating. 
Below  an  altitude  of  roughly  220  km,  the  electron  collision  frequency  increases  by  6%,  and 
is  self-similar  to  the  pump  shape.  Above  220  km,  the  electron  collision  frequency  changes 
due  to  temperatxire  modifications  by  a  decrement  of  1.5%. 

Figtire  8- 23a  displays  n«(s)  vs.  altitude  at  6  second  intervals.  The  density  modification 
normalized  to  the  equilibrium  density,  6ne(s),  is  plotted  vs.  altitude  at  6  second  intervals, 
and  is  shown  in  Figure  8-236.  Figure  8-23c  exhibits  6ne(s)  near  the  pump  center  as  a 
function  of  heating  time.  Above  an  altitude  of  220  km,  6n«  is  negative.  The  density  changes 
never  saturate,  in  agreement  with  Hansen^^^l,  who  heated  on  the  order  of  10  minutes,  and 
still  did  not  see  a  steady  state  achieved  while  heating  a  nighttime  ionosphere. 

The  behavior  of  6ne  is  better  tmderstood  upon  examination  of  Figures  8-24  and  8-25. 
Figiure  8-24  is  a  plot  of  the  recombination  rate  along  the  field  line,  /de(s),  vs.  altitude. 
Below  an  altitude  of  220  km,  there  is  a  5  -  10%  decrease  in  Above  an  altitude  of  220 
km,  there  is  roughly  a  10%  increase  in  Figtire  8-25  displays  the  diffusion  coefficient, 
D(s),  normalized  to  Above  an  altitude  of  160  km,  the  diffusion  coefficient  is  greater 
than  zero. 

For  this  nighttime  scenario,  at  an  altitude  of  220  km,  there  is  a  topological  change  in 
the  density  modification.  For  different  electron  density  and  temperature  equilibria,  this 
region  of  topological  change  will  move  to  a  different  altitude.  For  case  0704,  below  an 
altitude  of  220  km,  there  is  a  density  increase  primarily  caused  by  the  decrease  in  the 
recombination,  as  discussed  for  case  0112.  Since  in  the  nighttime  scenario,  the  piimp 
reaches  an  altitude  where  there  is  some  thermal  conductivity,  there  is  now  an  increase  in 
Te  at  altitudes  higher  than  the  pump.  This  catises  a  build-up  of  pressure  above  the  pump 
center,  and  thus  the  digging  of  a  thermal  density  cavity  as  discussed  by  Shoucri  et  al.^^^  to 
achieve  pressure  equilibrium. 

An  increase  in  n«  below  220  km  in  altitude,  and  a  decrease  in  ne  above  220  km  in 
altitude  means  that  at  220  km  in  altitude,  6n«  is  zero.  This  suggests  that  perhaps  the  best 
manner  in  which  the  radar  can  operate  is  to  align  the  null  point  in  6ne  near  the  caustic 
surface  by  tuning  the  antenna  to  the  appropriate  frequency.  This  should  minimize  the 
density  modifications  which  occur  in  the  ionosphere. 

This  null  point  in  6n«  remains  at  an  altitude  of  220  km  for  all  field  lines  on  all  meridizin 
planes.  When  the  pump  is  centered  at  an  altitude  significantly  lower  than  that  of  220 
km,  only  a  density  increase  is  observed  becatise  the  thermal  conductivity  is  vanishingly 
small.  As  the  HF  beam  moves  higher  in  altitude  ,  the  topological  alteration  in  the  density 
modification  is  observed. 

8-3.  Effect  of  Different  Recombination  Models 

The  transport  code  currently  has  the  capability  to  handle  three  different  recombination 
models.  The  first  is  that  which  has  already  beox  discussed,  namely,  that  the  recombination 
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rate,  /9e,  is  inversely  proportional  to  Te  raised  to  some  empirically  determined  power.  This 
model  causes  a  density  increase  when  the  pump  is  centered  on  a  field  line  at  a  low  altitude. 

The  second  model  is  that  of  Hansen^^^].  The  recombination  rate  remains  constant  in 
time,  but  is  a  function  of  position.  This  is  handled  in  the  transport  code  by  establishing  the 
recombination  rate  in  equilibrium,  but  never  updating  it  in  time.  When  this  recombinative 
model  is  utilized,  the  only  density  modifications  observed  at  all  altitudes  axe  thermal 
cavities  caused  by  pressure  imbalances.  In  the  lower  ionosphere,  such  density  cavities  are 
on  the  order  of  1%.  Figure  8-26  is  a  plot  of  ^Te(s)  vs.  alfirude  with  /Se(s,t)  =  /Se(s,t  =  0) 
for  case  0112.  Figure  8-27  displays  ^ne(s)  vs.  altitude  with  )3e(s,t)  =  /Se(s,t  =  0)  for  case 
0112. 

The  third  model  is  that  of  Shoucri  et  al.^^^;  i.e.,  there  is  not  a  recombinative  process 
explicitly  included  in  the  model.  The  recombination  term  jSne  is  lumped  in  with  the  source 
term  and  is  never  updated.  (In  the  model  of  Hansen^^^^,  n«  is  updated  in  the  term  /^eHe, 
but  fit  is  not.)  This  limit  is  achieved  by  setting  to  zero  within  the  code  the  molecular 
ion  densities,  as  then  fit  is  also  zero.  When  the  transport  for  this  limit  is  examined, 
it  is  found  that,  once  again,  density  modifications  at  all  altitudes  are  thermal  cavities. 
However,  the  cavity  does  not  saturate,  in  accord  with  Shoucri  et  al.^^l  Figures  8-28  and 
8-29  display  ^n«(8)  vs.  heating  time  when  the  recombination  rate  is  constant  in  time,  and 
when  the  recombination  rate  is  zero,  respectively.  Clearly,  when  there  is  recombination, 
the  modification  saturates. 

These  three  limits  are  realized  by  setting  a  switch,  “ibeta”,  to  the  appropriate  value 
in  the  input  file  “ht-oth.da”.  If  “ibeta”  =  0,  the  recombination  rate  is  updated  in  time;  if 
“ibeta”  =  1,  the  recombination  rate  is  constant  in  time;  if  “ibeta”  =  2,  the  recombination 
rate  is  zero. 
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Figure  8-2.  Unperturbed  Rays  in  the  Y-Z  Plane  for  Case  0112 


Figure  8-3.  Unperturbed  Rays  in  the  X-Y  Plane  for  Case  0112 
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Figure  8-4.  Unperturbed  Rays  in  the  X-Z  Plane  for  Case  0704 
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Figure  8<6.  Unperturbed  Rays  in  the  X-Y  Plane  for  Case  0704 
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Figure  8-7.  Typical  Q)|f  Profile  for  Case  0112 
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Figure  8-10.  Effect  of  a  Non-Sel^Coiuistent  Density 
Perturbation  of  20%  placed  in  the  E  repon  for  Case  0112 

(Y-Z  Plane) 


RAY  TRAJECTORY 


Figure  8-11.  Effect  of  a  Non-Self-Consutent  Density 
Perturbation  of  20%  placed  in  the  £  region  for  Case  0112 

(X-Y  Plane) 
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Figure  8-12.  Normalized  Thermal  Conductivity  vz.  Altitude 
for  Case  0112 
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Figure  8-13a.  Electron  Temperature  vs.  Altitude  for 
Case  0112 
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Figure  8-13b.  Normalised  Electron  Teiiq>erature  Modification 
va.  Altitude  for  Case  0112 
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Figure  8-13c.  Normalized  Electron  Teiiq>erature  Modification 
V8.  Time  Near  Peak  of  Q^f  Profile  for  Case  0112 
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Figure  8-14.  Electron  Cooling  vs.  Altitude  for  Case  0112 
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Figure  8-15.  Electron  Collision  Frequency  vs.  Altitude 
for  Case  0112 
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Figure  8>10a.  Electron  Density  vs.  Altitude  for  Case  0112 


ELECTRON  DENSITY 


62 


D*lta 


Figure  8-16b.  Normalized  Electron  Density  Modification  vs 
Altitude  for  Case  0112 
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Figure  8- 16c.  Normalized  Electron  Density  Modification  vs 
Time  Near  Peak  of  Qhf  Profile  for  Case  0112 
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Figure  8*17.  Recombination  Rate  vs.  Altitude  for  Case  0112 
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Figure  8>18.  Normalized  DiflFiizioii  Coefficient  vs.  Altitude 
for  Case  0112 
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Figure  8-19.  Normalized  Thermal  Conductivity  vs.  Altitude 
for  Case  0704  at  6  Second  Intervals  During  Heating 
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Figure  8-20a.  Electron  Temperature  vs.  Altitude  for 
Case  0704  at  6  Second  Intervals  During  Heating 
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Figure  8-20b.  Normalized  Electron  Temperature  Modification  vs. 
Altitude  for  Case  0704  at  6  Second  Intervals  During  Heating 


DELTA  TE 


1M  120  140  160  180  200  220  240  260  280  300  320  340 

Altitud*.  Ka 


69 


Dalit  Ta  Hatr  Patk  ya.  Tiaa 


Figure  8-20c.  Normalized  Electron  Temperature  Modification 
vs.  Time  Near  Peak  of  Q^r  Profile  for  Case  07u4 
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Figure  8-22.  Electron  Collision  Frequency  vs.  Altitude 
for  Case  0704  at  6  Second  Intervals  During  Heating 
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Figure  8*23a.  Electron  Density  vs.  Altitude  for  Case  0704 
at  6  Second  Intervals  During  Heating 
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Figure  8-23b.  Normalized  Electron  Density  Modification 
vs.  Altitude  for  Case  0704  at  6  Second  Intervals 
During  Heating 
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Figure  8-23g.  Normalized  Electron  Density  Modification 
vs.  Time  Near  Peak  of  Qhr  Profile  for  Case  0704 
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Figure  8-24.  Recombination  Rate  vs.  Altitude  for  Case  0704 
at  6  Second  Intervals  During  Heating 
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Figure  8-26.  Normalized  Electron  Temperature  Modification 
vs.  Altitude  for  Case  0112  at  6  Second  Intervals  During 
Heating  When  Recombination  Rate  Is  Constant  in  Time 
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Figure  8-27.  Normalized  Electron  Density  Modification 
vs.  Altitude  for  Case  0112  at  6  Second  Intervals  During 
Heating  When  Recombination  Rate  Is  Constant  in  Time 
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Figure  8-28.  Normalized  Electron  Density  Modification 
vs.  Time  Near  Peak  of  Qi,f  Profile 
for  Case  0112  When  Recombination  Rate  Is  Constant  in  Time 
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Figure  8-29.  Normalized  Electron  Density  Modification 
vs.  Time  Near  Peak  of  Q^f  Profile 
for  Case  0112  When  Recombination  Rate  Is  Zero 
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9.  FEEDBACK  OF  IONOSPHERIC  MODIFICATION  TO  THE  RAY 
TRACING  CODE 

The  transport  code  writes  to  an  output  file,  “eos_oth.ou”,  the  number  of  grid  points 
along  a  field  line,  and  the  size  of  the  grid  point  spacing;  the  ntimber  of  field  lines  in  a 
meridian  plane;  the  number  of  meridian  planes  on  which  the  modifications  in  density  and 
temperature  were  calculated;  the  changes  in  temperature,  ATe  =  Te(s,  t)  —  Te(s,  t  =  0); 
the  changes  in  density.  An*  =  ne(s,  t)  —  n€(s,  t  =  0);  the  beginning  point  of  each  field 
line,  So;  and  the  end  point  of  each  field  line,  Sf.  The  grid  point  spacing  and  the  nimiber  of 
grid  points  on  each  field  line  are  identical,  as  is  the  number  of  field  lines  in  each  meridian 
plane. 

The  ray  tracing  code  utilizes  this  information  in  the  following  manner.  At  each  point 
along  a  geometric  optics  ray,  an  evaluation  is  performed  to  determine  if  that  point  lies 
within  the  box  of  modification  specified  by  the  points  r^,  £2,  13,  and  £4,  which  were 
explained  in  Section  8.  If  the  point  is  outside  the  box  of  modification,  then  the  ray 
tracing  code  proceeds  as  usual.  If  the  point  lies  within  the  box  of  modification,  then  it 
is  necessary  to  find  the  value  of  ATe,  Ane,  the  first  derivatives  of  Ane,  and  the  second 
derivatives  of  Ane  &t  that  point.  A  numerical  interpolation  scheme  is  implemented  to  find 
these  quantities.  The  value  of  ATe  is  used  to  update  the  collision  frequency,  and  the  value 
of  Aue,  as  well  as  its  derivatives,  is  used  to  update  the  dielectric  tensor. 

The  interpolation  is  performed  in  the  magnetic  field  frame  of  reference.  Since  a  scalar 
quantity  does  not  depend  upon  the  reference  frame,  the  values  of  ATe  ^d  Ane  £ire  the 
same  in  the  laimch  point  frame  and  the  magnetic  field  frame.  However,  the  derivatives  of 
Ane  must  be  transformed  to  the  laimch  point  frame,  and  in  that  frame  they  are  given  by 
the  expressions 

frame  =  ^  *  (^^)b  frame 

(VVa)|,p.  S-'  ■  (VVn)g  ■  I  .  (2) 

Once  these  values  are  determined  and  the  dielectric  tensor  is  updated,  the  ray  tracing 
czdctdation  proceeds  as  before. 

9-1.  Da3rtiine  Heating  -  Case  0112 

It  is  worthwhile  to  examine  the  contour  plots  of  An^  and  its  derivatives  in  the  y-z 
plane  of  the  launch  point  reference  frame,  which  is  nearly  perpendicular  to  the  HF  beam. 
Figures  9-la-d  exhibit  the  behavior  of  these  quantities  in  the  box  of  modification. 

The  contours  of  Ane  m-e  positive,  and  centered  about  the  HF  beam  on  all  meridian 
planes.  The  derivatives  of  Ane  show  the  correct  topological  behavior,  namely  that  they 
are  positive  below  the  center  line  of  the  density  change,  and  negative  above. 

Figures  9-2  depict  the  paths  of  the  rays  after  heating  has  occurred.  The  introduction 
of  a  “blob”  of  density  to  the  ionosphere  acts  as  a  scatterer,  and  the  rays  deflect  outward 
both  in  elevation  and  azimuthal  angles.  The  point  of  entry  into  the  ionsphere  is  clearly 
observed  in  Figure  9-2a,  after  which  the  ray  paths  change  from  their  unperturbed  form. 
Figure  9-2c  shows  a  crossing  of  rays  in  the  x-y  plane.  The  rays  above  the  y=0  line  are 
deflected  upward,  and  those  below  the  y=0  line  ju-e  deflected  downward. 
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9-2.  Nighttime  Heating  -  Case  0704 

Figures  9-3  are  contour  plots  of  An®  and  its  derivatives  in  the  y-z  plane  of  the  launch 
point  reference  frame  over  the  range  of  the  box  of  modification.  At  lower  altitudes, 
the  density  modification  is  positive  definite  everywhere.  At  higher  altitudes,  the  density 
changes  are  negative  at  low  z,  and  positive  at  high  z,  as  discussed  in  Section  8. 

Figures  9-4  show  the  ray  paths  after  heating  has  taken  place  for  the  nighttime  scenario 
in  the  x-z,  y-z,  and  x-y  planes,  respectively.  The  ray  trajectories  in  the  x-z  plane  (Figure  9- 
4a)  chzmge  very  shghtly  from  their  vmperturbed  paths.  There  is  a  slight  spreading  present 
in  the  ray  paths  in  both  the  y-z  (Figure  9-45)  and  the  x-y  (Figure  9-4c)  plzines.  The  reason 
that  there  is  very  little  change  in  the  paths  of  the  rays  is  because  the  alterations  in  density 
and  temperature  which  occurred  are  very  small. 

9-3.  Iterative  Self-Consistent  Coupling 

While  this  analysis  only  goes  through  one  iteration,  all  of  the  tools  required  to  perform 
an  iterative  self-consistent  calciilation  have  been  developed.  The  procedure  outhned  in 
this  report  must  be  iterated  tmtil  convergence  in  the  ray  paths  is  achieved.  If  convergence 
csumot  be  attained,  then  the  duration  time  of  ionospheric  heating  (Athg^ting)  before 
feedback  is  performed  needs  to  be  shortened  until  convergence  occurs. 
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Figure  9>la.  Electron  Density  Modification  Contours 
in  Y'Z  Plane  for  Case  0112 
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Figure  9-lb.  Contours  of  X-Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0112 
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Figure  9-lc.  Contours  of  Y-Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0112 
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Figure  9- Id.  Contours  of  Z-Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0112 
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Figure  9-2b.  Ray  Paths  in  the  Y-Z  Plane  after  Heating 

for  Case  0112 
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Figure  9-2c.  Ray  Paths  in  the  X-Y  Plane  after  Heating 

for  Case  0112 
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Figure  9-3a.  Electron  Density  Modification  Contours 
in  the  Y-Z  Plane  for  Case  0704 
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Figure  0-3b.  Contours  of  X-Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0704 
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Figure  9>3c.  Contours  of  Y-Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0704 
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Figure  9*3d.  Contours  of  Z>Derivative  of  Electron 
Density  Modification  in  Y-Z  Plane  for  Case  0704 
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Figure  9-4c.  Ray  Paths  in  the  X-Y  Plane  after  Heating 

for  Case  0704 
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10.  PROPAGATION  EFFECTS  ON  RADAR  SIGNALS 

Ionospheric  heating  by  a  high  power  radar  Jbeam  may  result  in  strong,  time-dependent 
refraction  of  the  radar  beam  which  will  distort  the  radar  footprint  and  introduce  spurious 
doppler  shifts  in  the  radar  signals.  Self-focusing  of  the  radar  beam  in  the  ionosphere 
may  significantly  enhance  attenuation  of  the  radar  signal  through  ohmic  dissipation. 
Ducting  of  radiation  in  the  modified  ionosphere  resulting  from  localized  heating  may 
radically  alter  the  propagation  path.  Parametric  instabilities  may  also  have  important 
effects.  They  may  induce  ionospheric  radar  clutter  effects  due  to  backscattering  of  the 
primary  radar  signal  fiom  instability  product  waves,  may  enhance  phase  noise  in  the 
radar  signal  by  dielectric  fluctuation  effects  on  forward  propagation,  and  may  cavise 
additional,  anomalous  ionospheric  attenuation  of  the  primary  rzidar  signal.  Furthermore, 
the  parametric  instabilities  and  self-focusing  phenomena  may  interact  with  each  other  to 
further  enhance  their  respective  effects.  A  primary  goal  of  the  present  endeavors  is  to 
determine  the  conditions  under  which  such  nonlinear  phenomena  can  occur  in  ionospheric 
propagation,  with  special  attention  to  direct  ionospheric  heating  by  the  radar  signal,  and 
characterize  their  effects  on  radar  signals  and  system  performance. 

When  powerful  HF  wave  signals  are  obliquely  incident  on  the  ionosphere  two  important 
phenomena  can  occur:  transport  and  instabilities.  The  former  issue,  discussed  in  the 
previous  three  sections,  addresses  how  bulk  heating  (caused  by  high  incident  power) 
changes  the  profiles  of  the  backgroimd  plasma.  The  latter  pertains  to  the  formation  of 
plasma  structures  due  to  the  onset  of  instabilities  triggered  by  high  incident  power.  For 
HF  radzLT  parameters,  the  thermal  self-focusing  instability  is  the  most  likely  candidate. 

The  self-focusing  of  powerful  radio  waves  has  been  the  subject  of  numerous  analytical 
and  experimental  investigations.  During  vertical  HF  ionospheric  modification  experiments 
performed  by  Dimcan  and  Behnket^l  in  1978  and  by  others  later,  structtires  were  measured 
which  drift  through  the  illuminated  HF  beam  at  velocities  between  15  to  60  m/s  with 
a  periodicity  of  1  to  2  minutes  (i.e.,  with  a  transverse  scale  size  of  a  few  kilometers). 
According  to  Perkins  and  Goldman^^^,  in  the  F  layer  the  threshold  for  onset  of  the  self- 
focusing  instabihty  is  given  by 


where 


(1) 

(2) 


If  10  km  <  /  <  50  km,  then  for  Pc^p  =  3.48  x  10“^  W/m^,  a  value  typical  of  the 
amoimt  of  power  entering  the  F-layer  at  150  km  in  altitude  for  the  daytime  scenario 
(case  0112),  6.5  x  10~^  sec~^  <75  3.25  x  10“^  sec~^.  Thus  the  amoimt  of  time  for 
the  instability  to  develop  is  on  the  order  of  several  seconds  in  the  main  beam  region,  and 
is  shorter  at  higher  altitudes  where  the  electron  density  is  larger  and  ion-neutral  collision 
rates  are  snudler. 

For  the  nighttime  scenario,  (case  0704),  a  typical  power  flux  density  at  an  altitude 
of  150  km  near  the  c«iter  of  the  beam  is  Pc,F  =  ^•57  x  10“^  W /m^.  [Note:  While  the 
power  flux  density  of  the  nighttime  case  is  roughly  equivalent  to  that  of  the  daytime  case. 


98 


the  total  amoTint  of  power  crossing  the  meridian  plane  where  the  center  of  the  HF  beam 
is  at  150  km  in  altitude  is  660  kW  for  the  nighttime  case,  as  compared  to  462  kW  for 
the  daytime  case.  However,  the  nighttime  beam  is  slightly  wider  than  the  daytime  beam 
(4.5®  X  4.5®  as  compared  to  3®  x  3®),  because  the  frequency  of  the  nighttime  scenario  is  10 
MHz  as  compared  to  15  MHz  for  the  daytime  scenario,  but  the  size  of  the  antenna  is  the 
same  for  both  cases.  Thus  the  beam  area  is  larger  for  case  0704  than  for  case  0112,  which 
is  what  causes  the  power  flux  densities  of  the  two  cases  to  be  the  same  order  of  magnitude 
even  though  the  total  power  of  the  daytime  case  is  2/3  that  of  the  nighttime  case.]  If 
10  km  <  /  <  50  km,  then  8.0  x  10~®  sec“^  <  7  <  4.0  x  10~^  sec“^.  Thus  the  amoimt 
of  time  for  the  instability  to  develop  is  on  the  order  of  hours  in  the  main  beam  region. 
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11.  CONCLUSIONS 


11>1.  Developed  Capability 

The  resiilt  of  this  work  is  a  powerful  and  unique  capability  to  nvunerically  model  the 
oblique  propagation  and  heating  of  HF  waves  in  the  ionosphere.  Geometrical  optics  ray 
tracing  coupled  with  ionospheric  E  and  F  region  transport  models  can  self-consistently 
compute  the  nonlinear  deformation  of  the  HF  beam  due  to  heating  of  the  ionosphere  when 
the  waves  are  radiated  from  high  power,  high  gain  transmitters.  Such  treinsmitters  are 
fo\md  in  some  over-the-horizon  radar  systems  as  well  as  beamed  short  wave  broadcasting 
stations. 

The  work  reported  here  focused  on  developing  and  enhancing  capabilities  of  a 
three  dimensional  geometrical  optics  ray  tracing  code  that  can  propagate  rays  in  an 
inhomogeneus,  dissipative  and  anisotropic  medium  such  as  the  ionosphere.  The  code 
calculates  amplitude,  polarization  and  phase  of  the  wave  fields  along  the  rays,  so  that 
the  total  vector  field,  Poynting  flux  and  power  dissipated  at  a  given  point,  or  region,  in  the 
beam  can  be  calculated.  In  addition,  the  code  supports  time  dependence  of  the  medium. 
Doppler  shifts  caused  by  temporal  variations  in  the  medium,  such  as  drifts,  are  important 
in  measurements  of  radar  target  and  sea  surface  backscatter.  Such  drifts  can  be  modeled 
by  the  code. 

Similarly,  the  transport  code  enhancements  focused  on  realistic  modeling  of  the  physics 
of  ionospheric  heating  in  the  D,  E  and  F  regions.  In  order  to  be  useful  to  users  wanting  to 
model  a  specific  condition  of  the  ionsphere,  the  code  assumes  the  background  ionosphere 
to  be  a  given,  either  from  empirical  models  such  as  IRI,  measured  ionograms  at  the 
transmitter  site,  or  analytical  profiles  constructed  to  model  specific  phenomena  such  as 
drifting  striations.  The  code  then  calculates  the  deformation  of  the  backgrotmd  ionosphere 
caused  by  the  ohmic  dissipation  of  the  HF  beam.  Changes  in  electron  temperature  and 
density  modify  the  backgrotmd  equilibrium  which  in  turn  redistributes  the  ray  path. 
This  powerful  prescription  provides  the  most  complete  self-consistent  modeling  of  the 
propagation  and  heating  of  powerful  HF  waves  in  the  ionosphere  ctirrently  available. 

11-2.  Major  Findings 

In  order  to  demonstrate  the  capabilities  of  the  model,  sample  cases  were  chosen  to 
emphasize  specific  features  that  are  uniquely  handled  by  the  codes. 

Examples  of  ray  tracing  in  an  inhomogeneous  and  magnetized  plasma  were  shown. 
When  the  HF  beam  frequency  is  not  much  higher  than  the  electron  gyrofrequency  (i.e.,  5 
to  10  MHz)  and  the  beam  extends  well  into  the  region  of  high  plasma  frequency  (reflection 
close  to  f^F2),  the  effect  of  the  smisotropy  is  significant  and  the  beam  splits  into  two  distinct 
separate  modes  (sometimes  referred  to  as  ordinary  and  extraordinary  modes).  Under 
conditions  where  the  HF  beam  frequency  is  high  with  respect  to  the  electron  gyrofrequency 
or  the  plasma  frequency  is  small  compared  to  the  HF  beam  frequency  (reflection  in  the 
lower  ionosphere),  the  effect  of  the  anisotropy  is  negligible. 

Ducted  modes  are  also  of  interest.  Diurnal  variations  of  the  ionospheric  profiles  through 
different  longitudes  may  cause  a  beam  at  grazing  incidence  to  propagate  for  quite  long 
distances.  An  example  of  such  case  was  shown  in  our  results. 

A  daytime  heating  case  representative  of  a  mid-latitude,  high-simspot  number,  winter 
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conditions  was  tested  for  heating  using  a  powerful  transmitter  of  about  95  dBW.  The 
results  of  the  transport  analysis  show  that  about  60the  HF  energy  is  absorbed  ohmically 
in  the  lower  ionosphere  at  the  entrance  point  on  the  transmission  path.  Likewise,  when  the 
beam  comes  down  from  the  ionosphere,  high  losses  may  be  expected.  The  heating  of  the 
lower  ionosphere  causes  the  electron  temperature  to  increase  by  over  70%  zmd  the  density 
to  increase  by  about  20%.  Such  a  density  enhancement  can  cause  significant  additional 
refraction  and  defocusing  of  the  HF  beam. 

A  similar  summer  nighttime  case  shows  negligible  absorption  of  the  HF  energy  due 
to  greatly  reduced  collisionality  of  the  plasma..  Consequently,  electron  temperature  and 
density  changes  axe  also  negligible.  E-region  transport  is  mainly  local  due  to  high 
coUisionedity.  HF  heating  causes  an  increase  in  temperature  which  leads  to  an  increase 
in  density  due  to  reduced  recombination  rates.  In  the  F  region,  thermal  conduction 
tends  to  dominate  because  of  the  absence  of  molecular  ions  which  are  responsible  for 
recombination.  Increase  in  electron  temperature  induces  density  cavities  in  order  to  relax 
pressme  gradients. 

In  spite  of  heavy  daytime  absorption  in  the  E  region,  F-region  power  densities  are  in 
excess  of  theoretical  thresholds  for  thermal  self-focusing  instabilities.  For  the  nighttime 
case,  reduced  absorption  in  the  lower  ionosphere  is  offset  by  the  broader  beamwidth 
restdting  from  selection  of  lower  frequency  with  the  same  antenna  gain.  Therefore,  F- 
region  power  densities  in  the  nighttime  case  axe  compaxable  to  those  for  the  daytime  case. 

11-3.  Future  Work 

Continuation  of  the  work  reported  here  may  be  carried  out  in  a  several  directions. 
The  research  tools  developed  here  axe  useful  for  direct  simulation  of  HF  radar  systems. 
A  detailed  description  of  power  deposition  near  caustic  surfaces  and  in  strong  refraction 
zones  should  be  developed  to  transcend  the  hmitations  of  short-wavelength  asymptotics 
in  the  wave  propagation  model.  Some  straightforward  improvements  in  the  ray  tracing 
2dgorithm  can  lead  to  generally  coherent  phsise  information  associated  with  radeir  waves, 
and  this  wo\ild  allow  the  direct  simulation  of  ionospheric  effects  on  target  and  sea  clutter 
doppler  spectra  in  HF  radars.  Studies  of  propagation  through  striations  may  be  carried 
out  for  realistic  striation  patterns.  Cross-field  tramsport  caused  by  drifting  striations  can 
be  accoimted  for  by  adding  two  dimensional  effects  to  the  transport  analysis.  Finally,  a 
systematic  numerical  study  of  the  thermal  self-foc\ising  instability  could  be  carried  out  to 
explore  and  validate  existing  anedytical  results  on  the  subject. 
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APPENDIX  A.  COLLISIONAL  RATES 

The  following  is  a  svimmary  of  the  coUisional  rates  used  in  the  transport  model  with 
reference  to  where  they  were  obtained. 

1.  Electron-Neutral  Collision ; 

a.  (Ref.  23,  25) 

i/e,Nj  =  2.33  X  10~^*n(N2)Il  -  1.21  x  10"^Te]Te  s'^  . 

b.  (Ref.  23,  25) 

i/e,Oj  =  1.82  X  10”^VO2)[l  +  .036Ty^]Ty^  . 

c.  (Ref.  24) 

i/e,o  =  2.65  X  10-^VO)Ty^  s-1  . 

d.  (Ref.  23,  25) 

i/e,He  =  4.6  X  10“^  VHe)Ty^  s'^  . 

e.  (Ref.  23,  25) 

i/e,H  =  4.35  X  10" VH)[1  -  1.35  X  lO'^^TelXy^  s"^  . 


2.  Ion-Neutral  Collisions 

a.  (Ref.  23,  25) 

1/0+ Nj  =  6.82  X  10-^0n(N2)  s'^  . 

b.  (Ref.  23,  25) 

*'0+,02  =  ^  10"^®n(O2)  . 


c.  (Ref.  23,  25) 

^0+,0  =  3.67xlO-lln(0)[(Ti-|-Tn)/2]l/2^1.0-.0641ogio[(Ti-hTn)/2]}2  s"^  . 


=  1.32  X  10"^VHe)  s"l  . 


d.  (Ref.  23,  25) 


*'0+,He 


e.  (Ref.  25) 


*^0+3  =  4.13  X  10-1  VH)Ty^[1.0  -  .0471ogio(Ti)]2  s"!  . 


f.  (Ref.  23,  25) 

=  4.13  X  IO-IVN2)  s-1  . 

g.  (Ref.  25) 

''OJ.O,  = 

2.59  X  10-lln(O2)[(Ti  -f  Tn)/2]l/2 

{1.0-.0731ogio[(Ti-hTn)/2]}2  s-1  . 

h.  (Ref.  23,  25) 

• 

=  2.31  X  10-l®n(O)  s-1  . 

i.  (Ref.  23,  25) 

=  7.0  X  lO-lln(He)  s"!  . 

j.  (Ref.  23,  25) 

"'Ot.H 

=  6.5  X  lO-lln(H)  s-1  . 

k.  (Ref.  23,  25) 

‘'NO+.Nj 

*  4.34  X  IO-IVN2)  s-1  . 

1.  (Ref.  25) 

*^0+.0} 

=  4.27  X  10-l°n(O2)  s"!  . 

m.  (Ref.  23,  25) 

''NO+.O 

*  2.44  X  10-l°n(O)  s"!  . 

n.  (Ref.  23,  25) 

*^0+.Hp 

=  7.4  X  10-1  In(He)  s"!  . 

o.  (Ref.  23,  25) 

"NO+.H 

=  6.9  X  lO-lln(H)  s-1  . 

3.  Electron-Ion  Collisions 


a.  (Ref.  23,  24,  25) 


54.5ni 

t3/2 
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APPENDIX  B.  ELECTRON  COOLING  RATES 

The  following  is  a  summary  of  the  electron  cooling  rates  used  in  the  transport  model 
with  reference  to  where  they  were  obtained. 

A.  Cooling  due  to  translational  electron-neutral  interactions 

1.  (Ref.  24,  25) 

Le(e,N2)  =  1.77  xlO“^^nen(N2)[l- 1.21  xlO“'*Te]Te(Te-Tn)  . 

2.  (Ref.  24,  25) 

Le(e,02)  =  1.21  X  10-^®nen(O2)(l  +  .036Ty^]Ty^(Te  -  Tn)  eVcm-^3-^  . 

3.  (Ref.  25) 

Le(e,0)  =  7.9  X  10-^®nen(O)Ty^(Te  -  Tn)[l -H  5.7  X  10-'‘Te]  eVcm~^s-^  . 

4.  (Ref.  24,  25) 

Le(e,H)  =  9.63  xlO“^®nen(H)Ty^(Te-Tn)[l- 1.35  xlO-^Te]  tVcrrC^s-^  . 

5.  (Ref.  24,  25) 

Le(e,He)  =  2.46  x  10**“nen(He)Ty^(Te  -  Tn)  eVcm-^s-^  . 

B.  Cooling  due  to  rotational  excitatior  of  N2  and  O2 

1.  (Ref.  24,  25) 

Le(e,N2rot)  =  2.9  x  10-^Vn(N2)(Te  -  Tn)Te"^/^  eVcm'h'^  . 

2.  (Ref.  23,  24,  25) 

Le(e,02rot)  =  6.9  x  10-^'‘nen(O2)(Te  -  Tn)Tr^''^  eVcm-^^-l  . 

C.  Cooling  due  to  vibrational  excitation  of  N2  and  O2 

1.  (Ref.  25) 

Le(e,N2vib)  =  2.99  x  10"^2n^n(N2)exp[f(Te  -  2000)/(2000Te)] 
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X {[1  -  exp[g(Tn  -  Te)/(TeTn)]}  eVcm-h-^ 
f  =  1.06  X  10'*  +  7.51  X  102tanh[l.l  x  10"3(Te  -  1800)] 
g  =  3300  +  1.233(Te  -  1000)  -  2.056  x  10“'*(Te  -  1000)(Te  -  4000). 

2.  (Ref.  25) 

Le(e,02vib)  =  5.196  x  10-*2n«n(O2)exp[h(Te  -  700)/(700Te)] 

x{l-exp[2770(Tn-Te)/(TeTn)]}  eVcm-^3-^ 
h  =  3300  -  839sin[1.91  x  10~'*(Te  -  2700)]. 

D.  Cooling  due  to  fine  struct  tire  excitation  of  0 

1.  (Ref.  25) 

Le(e,O^P)  =  8.629  x  Ai(B!)iTef®‘"*/2l{ei(Di  -  Ei)  +  5.91  x  10" 

^  f=l 

x(T„  -  Te)[(l  +  Bi)Di  +  (^^  +  1  +  Bi)Ei]}  cFcm-^s-i 

■te 

where 

€i  =  .02(i  =  l),  .028  (i- 2),  .008  (i  =  3) 

Di  =  exp[-228/Tn]  (i  =  1),  exp{-326Tn]  (i  =  2),  exp[-326Tn]  (i  =  3) 

Ei  =  exp[-228/Te]  (i  =  1),  exp[-326Te]  (i  =  2), 
exp[-98/Te  -  228Tn]  (i  =  3) 

(EE)i  =  228  (i  =  1),  326  (i  =  2),  98  (i  =  3) 

Ai  =  8.58  X  10"®  (i  =  1),  7.201  x  10"®  (i  =  2),  2.463  x  10"^  (i  =  3) 

Bi  =  1.019  (i  =  1),  .8998  (i  =  2),  1.268  (i  =  3) 

(B!)i  =  1.019  (i  =  1),  .9169  (i  =  2),  1.1626  (i  =  3) 

with 

Z  =  5  +  3e"22®/T„  +  e"326/T„  . 

E.  Cooling  due  to  electronic  excitation  of  O  and  O2 

1.  (Ref.  25) 

Le(e,  0(*D))  =  1.57  x  10"*2n^n(O)exp[d(Te  -  3000)/(3000Te)] 
x{l  -exp[22,713(Tn  ~Te)/(TeTn)]}  cVcm"3s"*  , 
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where 


d  =  2.4  X  lO"**  +  0.3(Te  -  1600)  . 


2.  (Ref.  25) 

Le(e,02(^Ag))  =  1.143  x  10“^^nen(O2)exp[f(Te  -  1500)/(1500Te)] 
x{l  -exp[ll,400(Tn-Te)/(TeT„)I}  eVcm'h-^  , 

where 

f  s  {13,200  +  1410sin[.000241(Te  -  500)}  x  {1  +  exp[(Te  -  14,011)/1048]} 

3.  (Ref.  25) 

L,(e,Oj(‘Eg))  =  1.616  x  10-*®n,n(O2)e)cp[g(Te  -  1600)/(1500T,)1 
x{l  -otp(ll,400(T„  -T.)/(TeT„)l}  eVcm-’j'*  , 

where 

g  s  {19,225  +  560sm[.000383(Te  -  1000)}  x  {1  +  exp[(Te  -  16,382)/1760]} 
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APPENDIX  C.  CONFERENCE  PRESENTATIONS 

The  following  presentations  and  publications  have  resiilted  from  contract  sponsorship. 

1.  Modeling  of  High  Frequency  Propagation  and  Heating  in  the  Ionosphere,  AGARD 
Conference  on  Ionospheric  Modification  and  Its  Potential  to  Enhance  or  Degrade 
the  Performance  of  Mihtary  Systems,  Electromagnetic  Wave  Propagation  Panel 
Symposium,  Bergen,  Norway,  28  -  31  May,  1990. 

2.  Modeling  of  High  Frequency  Propagation  and  Heating  in  the  Ionosphere,  in  AGARD 
Conference  Proceedings  No.  485  on  Ionospheric  Modification  and  Its  Potential  to 
Enhance  or  Degrade  the  Performance  of  Military  Systems,  NTIS,  Springfield,  VA, 
1990. 

3.  Oblique  Ionospheric  Heating  by  HF  Waves,  1991  North  American  Radio  Science 
Meeting,  Commision  H  -  Ionospheric  Heating  and  Modification,  London,  ONT, 
Canada,  J\me  24  -  28,  1991. 

4.  Oblique  Ionospheric  Heating  by  HF  Wai^s,  Annual  Meeting  of  the  American  Physical 
Society,  Division  of  Plasma  Physics,  Tampa,  FL,  November  4-8,  1991. 

5.  Modeling  of  High  Frequency  Propagation  and  Heating  in  the  Ionosphere,  Over- 
the-Horizon  Radar  Technology  Conference,  Electronic  Systems  Division  and  Rome 
Laboratories,  Hanscom  AFB,  MA,  November  5-7,  1991. 
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APPENDIX  D.  COMPUTER  CODES 


To  obtain  further  information  concerning  the  computer  codes  (Appendix  D) 
for  this  report,  contact: 

RL  (ERCP/ Stanford  Yukon) 

Hanscom  AFB  MA  01731-5000 

(617)  893-5353 
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MISSION 


OP 

ROME  LABORATORY 

Rome  Laboratory  plans  and  executes  an  interdisciplinary  program  in  re¬ 
search,  development^  test,  and  technology  transition  in  support  of  Air 

3 

Force  Command^  Cwitrol,  Communications  and  Intelligence  (C  I)  activities 
for  all  Air  Force  platforms.  It  also  executes  selected  acquisition  programs 
in  several  areas  of  expertise.  Technical  and  engineering  support  within 
areas  of  competence  is  provided  to  ESD  Program  Offices  (POs)  and  other 
ESD  elements  to  perform  effective  acquisition  of  C  l  systems.  In  additiont 
Rome  Laboratory's  technology  supports  other  AFSC  Product  DivisioMf  the 
Air  Force  user  community,  and  other  DOD  and  nonrDOD  agencies.  Rome 
Laboratory  maintains  technical  competence  and  research  programs  in  areas 
including,  but  not  limited  to,  communications,  command  and  control,  battle 
management,  intelligence  informatiwi  processing,  computational  sciences 
and  software  producibility,  wide  area  surveillance/sensors,  signal  proces¬ 
sing,  solid  state  sciences,  photonics,  electromagnetic  technology,  super¬ 
conductivity,  and  electronic  reliability/maintainability  and  testability. 


